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20.  ABSTRACT  (Continued) 

$ection(C)  Fluid  Dynamics  and  Ablative  Acceleration' 

Secti^(D)  Development  of  Computational  Techniques^ 

^ctioh(E^  Rayleigh-Taylor  Stabilization  Techniques  , 

Our  work  on  the  transport  of  charged  particle  beams  has  concentrated  on  two  problems;  (Al)  deter- 
mining computationally  the  emittance  growth  due  to  plasma  effects  from  space  charge  forces  in  heavy 
ion  beam  transport  through  a quadrupole  focussing  lens  system  and  (A2)  ion  beam  transport  and  sta- 
bility in  preprepared  plasma  channels.  In  Section  B we  report  our  work  on  the  absorption  of  laser 
light  in  the  plasma  blow  off  at  material  surfaces.  The  role  of  non-resonant  absorption  due  to  ion 
acoustic  fluctuations  has  been  studied  in  the  presence  of  magnetic  fields.  The  thrust  of  our  work  on 
fluid  dynamics,  laser  ablative  acceleration,  ion  deposition  and  computational  techniques  is  focussed  on 
the  scoping  and  subsequent  quantitative  analysis  of  Rayleigh-Taylor  instability  problems  and  closely 
allied  asymmetry  problems.  The  occurrence  or  potential  occurrence  of  growing  R-T  modes  in  pellet 
implosion  is  nearly  universal.  Nevertheless,  four  rather  distinct  environments  can  be  identified  where 
the  operative  physics  and  hence  the  required  computational  tools  benefit  from  independent  analysis. 
The  organization  of  the  remaining  four  sections  of  the  annual  report  highlights  these  distinct  environ- 
ments. Each  of  these  four  segments  combines  elements  of  C,  D and  E above. 

C.  Ablation  Layer  Structure  and  Analysis  of  Rayleigh-Taylor  Instability. 

D.  Rayleigh-Taylor  Instability  at  Internal  Material  Interfaces. 

E.  Instability  and  Symmetry  Issues  in  Colliding  Shells. 

F.  Rayleigh-Taylor  Stability  and  Central  Core  Dynamics. 

Each  of  these  four  aspects  differs  enough  from  the  others  to  make  treating  them  separately  worthwhile. 
The  Rayleigh-Taylor  Instability  is  an  ideal  unifying  theme  for  much  of  the  research  which  has  to  be 
done. 

The  report  concludes  with  a small  section  on  our  plans  for  the  extension  and  continuation  of 
this  research  and  an  NRL  Inertial  Confinement  Fusion  bibliography  encompassing  recent  contributions 
to  ICF  theory  and  understanding.  Footnotes,  figures,  and  references  for  each  section  below  are  com- 
piled separately. 
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Al.  Space  Charge  Limited  Transport  of  Heavy  Ion  Beams 

Current  proposals  for  heavy  in  beam  ignited  inertial  fusion  are  attractive  because  they 
generally  require  only  modest  extensions  of  current  technology.  However,  the  ability  to  focus 
the  required  power  onto  a small  spot  at  a distance  from  the  final  lens  requires  a beam  not  only 
with  high  power  but  also  low  emittance.  It  has  been  known  for  some  time  from  the  observed 
behavior  in  proton  linacs  and  from  the  theoretical  analysis  of  beam  stability  that  high  phase 
space  current  densities  are  accompanied  by  space-charge  driven  instabilities  and  possible  emit- 
tance growth.  Because  emittance  growth  in  turn  limits  the  spot  size  to  which  a beam  can  be 
focused,  these  instabilities  can  create  a m^or  barrier  to  the  ignition  of  a pellet. 

When  the  current  in  a particle  beam  is  pushed  up  towards  its  space-charge  limit,  collective 
oscillations  that  can  be  supported  by  the  beam  are  plasma-like  in  behavior.  This  means  that  all 
the  non-linearities  and  attendant  complexities  of  the  Vlasov  equation  can  become  important  and 
render  a detailed  theoretical  treatment  extremely  complex.  Because  the  beam  is  only  a few 
Debye  lengths  across  and  is  transported  for  only  a few  plasma  periods  the  detailed  behavior  is 
particularly  amenable  to  computer  simulation  by  the  particle  techniques  well  developed  in  the 
plasma  physics  community.  Since  the  beams  will  be  long  compared  to  their  transverse  dimen- 
sions, the  transverse  behavior  is  simulated  in  two  dimensions  in  a reference  frame  moving  with 
the  beam.  The  orbits  of  several  thousand  rod-shaped  particles  are  followed  as  these  particles 
move  in  the  combination  of  lens  focusing  forces  and  their  self-consistently  calculated  space 
charge  forces. 
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Initial  calculations  have  been  done  assuming  a thin  lens,  alternate  gradient  transport  sys- 
tem. In  order  to  make  contact  with  existing  theory  and  the  ongoing  program  at  LBL,  an  initial 
Kapchinskii-Vladimirskij  (K-V)  distribution  was  assumed  for  the  particles.  The  lens  force  are 
simulated  by  applying  to  each  particle  a IkVx  and  A Vy  proportional  to  the  lens  strength  and  to 
the  X and  y displacement  of  each  particle  from  the  beam  axis.  In  the  absence  of  space  charge 
this  numerical  system  reproduces  the  idealized  thin  lens  system  to  within  roundoff. 


; To  ensure  accurate  and  reliable  simulations,  an  extensive  series  of  numerical  tests  were 

performed.  A test  case  was  chosen  to  be  an  alternate  gradient  transport  system  with  90"  phase 
advance  depressed  by  space  charge  to  30°.  A numerical  regime  was  identified  where  the  results 

- i 

' are  insensitive  to  variations  over  factors  of  two  in  time  step,  numbers  of  simulation  particles, 

number  of  computational  grids,  boundary  condition  and  macroparticle  size'.  By  extending 
these  results,  and  by  further  testing,  it  was  determined  that  simulation  could  accurately  be  per- 
formed down  to  a time  depression  of  at  least  (7.5°.) 

In  addition,  a comprehensive  series  of  tests^  was  performed  to  examine  the  long  term 
growth  in  emittance  due  to  the  smaller  than  physical  number  of  particles  in  the  simulation  sys- 
tem. A mstior  effect  was  seen  only  when  the  number  of  particles  was  decreased  to  far  below  the 
number  which  could  be  economically  employed. 

With  the  reliability  of  the  numerical  techniques  established,  a detailed  parametric  study 
was  undertaken  of  a 90°  thin  lens  quadrupole  transport  system.  Good  qualitative  agreement 
was  obtained  with  theoretical  calculations  at  LBL.  One  case^,  with  the  tune  depressed  to  45°, 
was  subject  to  detailed  analysis  and  quantitative  agreement  obtained  with  the  LBL  theory  with 
respect  to  growth  rates  and  the  detailed  form  of  the  unstable  behavior. 

Figure  Al.l  is  a plot  of  the  xand  y rms  emittances  of  a beam,  with  enough  space  charge 
to  depresr  the  tune  to  30°,  as  it  evolves  from  the  initial  K-V  equilibrium.  The  behavior  is 
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characterized  by  an  initial  constancy  of  the  emittance  followed  by  a period  of  rapid  growth.  The 
growth  rapidly  saturates  and  the  emittance  is  then  constant  for  the  duration  of  the  run.  This 
saturation  and  steady  state  is  also  observed  in  other  moments  of  the  distribution  function  and 
the  system  does  indeed  appear  to  have  reached  a steady  state.  For  any  given  equilibrium  distri- 
bution, Maschke  and  Courant^  have  derived  a scaling  law  for  transportable  power 

14/3 

P-C  j (1) 


j By  using  the  figure  of  merit  defined  by  Laslett  and  degrading  it  for  emittance  growth  and  beam 

I * 

i I spreading,  it  can  be  argued^  that  the  case  shown  in  Fig.  1 corresponds  to  C - 1.2  x 10'^  . 

[ '! 

^ . This  steady  state,  therefore  turns  out  to  yield  a transportable  power  consistent  with  pellet  igni- 

; i 

I ' * tion  requirements. 

When  the  current  being  transported  results  in  a much  greater  tune  depression,  the 

: behavior  of  the  emittance  evolution  changes.  Figure  A 1.2  shows  the  evolution  in  emittance  for 

r 

I , a tune  of  90°  depressed  by  space  charge  to  7.5°.  The  period  of  constancy  and  rapid  growth  in 

this  case  is  followed  by  a period  of  slower  growth  in  the  emittance.  A case  with  the  tune 
depressed  to  15°  shows  similar  behavior  but  the  long  term  growth  is  somewhat  slower.  These 
results  suggest  that  a space  charge  limit  does  exist  for  quadrupoie  transport,  and  when  it  is 
exceeded  the  emittance  will  grow  until  the  phase  space  density  has  been  reduced  to  below  this 
limiting  value. 

How  this  limiting  value  varies  as  a function  of  various  transport  system  parameters  has 
not  yet  been  established.  However,  a simulation  of  a 130°  transport  system  depressed  to  50* 
verified  the  predicted  instability  in  the  envelope  equations  and  showed  a rapid  growth  in  emit- 
tance by  about  a factor  of  five.  This  lends  support  to  the  assertion^  that  transport  systems  with 
advances  greater  than  90°  are  to  be  avoided.  In  fact  recent  work  by  Lasslett  and  Smith*  indi- 
* , cate  that  60°  transport  systems  hold  the  promise  of  transporting  even  greater  currents  than  the 

90°  systems.  Preliminary  simulations  appear  to  support  this  conclusion  at  least  for  K-V  distri- 
butions. 
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In  order  to  determine  the  dependence  of  space-charge  limited  transport  on  the  form  of 
the  initial  distribution  function,  a series  of  simulations  have  been  performed  with  different  ini- 
tial distributions.  In  order  to  overcome  the  difficulty  of  establishing  self  consistent  equilibria  in 
the  presence  of  space  charge,  these  distributions  are  established  at  low  currents.  The  currents 
are  then  slowly  increased  so  that  the  distribution  can  adjust  itself  in  a manner  similar  to  what 
would  occur  in  the  center  of  a beam  being  slowly  bunched  longitudinally.  The  form  of  the  dis- 
tribution does  in  fact  effect  both  the  growth  rates  and  instability  thresholds.  In  general 
however,  behavior  worse  than  the  K-V  results  has  not  been  observed  in  these  simulations. 


Progress  thus  far  has  shown  that  particle  simulations  are  a very  promising  approach  to 
reliably  modelling  intense  beam  behavior.  The  results  summarized  here  have  indicated  that 
space  charge  collective  effects  can  have  a significant  effect  on  beam  transport  and  focusing. 
Since  system  economics,  performance  and  even  feasibility  can  be  influenced  by  these  effects, 
and  the  numerical  model  employed  can  easily  incorporate  many  of  the  complications  important 
to  the  performance  of  an  actual  system,  computer  simulation  appears  to  be  a highly  cost 
effective  approach  in  making  intelligent  decisions  in  the  design  of  future  experiments. 
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A2.  Stability,  Channelling,  and  Transport  of  Ion  Beams 


With  recent  advances  in  the  technology  of  intense  light-ion  beam  generators,'-^  the 
suggestion  tha.  these  ion  beams  could  be  used  in  a pellet  fusion  scheme  has  drawn  considerable 
interest^  as  have  similar  concepts  using  heavy  ion  beams^  (e.g.  20  GeV  U*''  ions).  In  a func- 
tioning reactor,  the  ion  beam  would  have  to  be  propagated  a distance  L (about  10  m)  from  the 
acceleration  region  while  employing  some  focussing  scheme  to  deliver  the  beam  to  the  target. 
An  investigation  of  the  equilibrium  and  stability  of  such  ion  beams  is  presently  underway  in 
order  to  determine  the  feasibility  of  propagating  and  focussing  them  onto  a pellet. 

For  heavy  ion  beams,  field-free  ballistic  focussing  to  the  target  has  been  proposed.*-^  For 
light  ion  beams  (e.g.  10  MeV  protons),  one  concept^  that  is  under  investigation  involves 
focussing  the  beam  just  outside  the  diode  region  by  known  techniques'  and  then  iiyecting  the 
focussed  beam  into  a z-pinch  plasma  channel.  The  beam  remains  focussed  during  propagation 
by  the  magnetic  field  in  the  plasma  channel  and  is  delivered  on  target  with  no  additional  focuss- 
ing. For  both  of  these  techniques,  it  must  be  demonstrated  that  the  self-consistent  electromag- 
netic fields  associated  with  the  plasma  response  to  beam  passage  do  not  disrupt  transport  and 
focussing.  We  have  started  to  investigate  this  problem  by  concentrating  on  the  physics  occur- 
ring in  the  plasma  channel  as  a result  of  the  ion  beam  passage.  This  area  of  investigation  has 
been  largely  ignored  in  studies  so  far;  it  has  been  generally  assumed  that  the  channel  does  not 
respond  to  the  beam  on  the  time  scale  involved.  With  beam  pulse  lengths  of  SO  ns  for  light  ion 
beams  and  of  around  10  ns  or  more  for  heavy  ion  beams,  this  assumption  is  not  generally 
correct  for  transport  in  channels  of  cm  radius  or  less.  Magnetic  nonneutralization  effects  are 
important  for  large  radius  beams  as  well,  since  small  (i.e.  ~ 1 kA)  net  currents  are  sufficient  to 
destroy  ballistic  focussibility  to  mm  spot  sizes  for  heavy  ion  beams. 

Recent  work  has  concentrated  on  the  plasma  response  to  light-ion  beams  with  a view  to 
answering  the  questions: 
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(a)  Does  a plasma-parameter  window  exist  which  allows  efficient  transport  over  10  meter 
distances? 

(b)  Do  the  self-consistent  electromagnetic  fields  associated  with  beam  heating,  net 
current  effects  and  plasma  hydrodynamics  disrupt  beam  transport  or  its  ability  to  bunch? 

Some  constraints  which  have  to  be  satisf  ed  by  the  beam  plasma  system  include: 

(a)  The  density  in  the  hole  bored  by  the  applied  z-discharge  and  beam  currents  must  be 
sufficiently  high  to  permit  complete  charge-neutralization. 

(b)  The  beam-heated  plasma  must  be  sufficiently  high  in  temperature  to  provide  good 
magnetic  neutralization. 

(c)  The  temperature  must  be  sufficiently  high  to  reduce  return-current  electric  fields  to  a 
level  which  neither  slows  the  beam  appreciably  nor  disrupts  bunching. 

A 1-D  radial,  time-dependent,  Eulerian  code  has  been  developed  recently  to  answer  these 
questions.  The  code  solves  the  plasma  moment  equations.  Maxwell  equations,  and  transport 
equations  using  coefficients  for  electrical  and  thermal  conduction  taken  from  Braginskii.^  The 
beam  is  specified  by  its  current  density  in  space  and  time.  Direct  beam  energy  deposition  and 
radiation  cooling  as  well  as  return-current  heating  have  been  included  in  the  code.  Initia<  con- 
ditions consist  of  a plasma  z-discharge  in  a high-density  gas^.  The  density  of  the  discharge 
channel  is  typically  .01-.  1 atmosphere  of  low-z  gas  while  the  background  density  may  be  much 
higher. 

Preliminary  results  appear  in  Fig.  A2.1.  They  correspond  to  a 1 MA,  10  MeV  proton 
beam,  0.6  cm  in  radius  propagating  in  a preformed  hydrogen  channel  of  the  same  dimension 
with  an  initial  density  on  the  axis  of  2 x 10 Results  shown  a.'e  radial  profiles  for  density. 
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azimuthal-magnetic  and  axial-electric  fields,  and  temperature  at  r — 40  ns  into  the  beam  pulse. 
These  are  the  conditions  which  develop  near  the  tail  of  the  beam.  No  impurities  have  been 


!■ 
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included  in  this  calculation  so  that  only  hydrogen  bremsstrahlung  losses  are  included.  Time 
histories  through  any  point  fixed  in  the  chamber  of  the  maximum  B-field,  maximum  tempera- 
ture, central  electric-field  strength  and  minimum  channel  density  appear  in  Fig.  A2.2. 

Conclusions  which  can  be  drawn  at  this  point  are: 

(1)  Hydrogen  backgrounds  of  1-3  x 10'^;m~^  in  the  channel  allow  propagation  satisfying 
all  constraints  to  transport  and  bunching. 

(2)  Hydrodynamic  expansion  at  lower  densities  may  result  in  insufficient  neutralizing 
densities. 

(3)  Low  plasma  temperatures  at  higher  densities  cause  excessive  return-current  electric 
field  strengths  which  slow  the  beam  and  perturb  bunching.  Line-radiation  cooling  from  impuri- 
ties ih  the  backgiound  may  result  in  too  low  temperatures  if  present  in  sufficient  densities. 

Phenomena  which  have  not  been  included  in  the  code  yet  and  which  are  important  to  a 
more  complete  understanding  of  the  channel  physics  are: 

• Non-equilibrium  ionization  of  the  plasma  at  the  head  of  the  beam. 

• Values  for  electrical  conductivity  at  low  temperatures  which  determines  the  penetration  of 
the  B-held  into  the  plasma  and  the  degree  of  current  non-neutralization  achieved  late  on 
in  the  pulse. 

These  phenomena  plus  comparison  with  experiments  at  lower  current  intensities  (with 
correspondingly  lower  temperatures  and  ionization  levels),  extension  of  this  work  to  heavy  ion 


10 


beams,  and  a more  consistent  treatment  of  the  ion-current  distribution  in  position  and  time  due 
to  seif  fields  and  bunching  will  constitute  future  developments  of  this  work. 

During  the  past  year,  an  analysis  of  microinstabilities  associated  with  light  and  heavy  ion 
beams  traversing  the  plasma-background  conditions  established  by  the  MHD  code  was  carried 
out. 

The  linear  growth  of  two  electrostatic  microinstabilities,  the  e-b  mode  and  the  e-i  mode, 
have  been  investigated  under  the  conditions  of  good  charge  and  current  neutralization  for  both 
heavy  ion  beam-plasma  systems’  and  light  ion  beam-plasma  systems.  ” The  e-b  mode  is  driven 
by  the  relative  streaming  between  beam  ions  and  electrons  and  the  e-i  mode  is  driven  by 
streaming  between  plasma  ions  and  electrons.  If  either  of  these  two  modes  is  sufficiently 
unstable  to  create  strong  microturbulence,  the  resulting  anomalous  resistivity  could  seriously 
affect  the  current  neutralization  of  the  beam  which,  in  turn,  might  disrupt  beam  transport. 

For  heavy  ion  beams,’  it  was  found  that  both  modes  can  be  stabilized  in  parameter 
regimes  which  are  consistent  with  the  conditions  for  charge  and  current  neutrality  and  which 
are  appropriate  for  a fusion  reactor.  Typical  marginal  stability  curves  for  both  the  e-b  mode  and 
the  e-i  mode  are  shown  in  Fig.  A2.3  for  a 40  tieV  beam  of  U*'  ions  and  in  Fig.  A2.4  for  a 7 
GeV  beam  of  ions.  If  occurring  at  all,  these  instabilities  can  grow  only  when  the  beam  is 
quite  close  to  the  pellet.  Thus,  it  is  unlikely  that  either  of  these  electrostatic  modes  will 
significantly  degrade  beam  transport  or  focussing. 

For  light  ion  beams,^  the  high  densities  required  for  the  background  z-pinch  plasma  pro- 
vide for  a high  frequency  of  collisions  between  the  electrons  and  plasma  ions.  This  results  in 
the  collisional  damping  of  the  e-b  mode  for  T,  < 50  eV  in  a typical  system  with  the  z-pinch 
current  flowing  on  the  surface  of  the  plasma  channel.  If  the  z-pinch  current  flows  in  the  inte- 
rior of  the  plasma  channel,  the  azimuthal  magnetic  held  further  stabilizes  the  e-b  mode  < 


1 keV  for  stability).  Since  it  is  expected  that  ihe  current  distribution  will  be  between  these 
extreme  forms,  and  hydrodynamic  computer-code  calculation^  demonstrate  that  the  electrons 
will  not  be  heated  above  SO  eV,  the  e-b  mode  will  be  stable.  The  e-i  mode  is  also  stabilized  by 
collisional  damping.  Other  generally  more  slowly  growing  modes,  such  as  electromagnetic 
pinching  instabilities,  are  presently  under  investigation. 
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Figure  A2.1.  Radial  profiles  of  density,  azimuthal  magnetic  field 
cucial  electric  field,  and  temperature  at  t = 40  ns. 
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Figure  A2.3  Typical  marginal  stability  curves  as  functions  of  n^/i^,  the 
ratio  of  beam  radius  to  initial  beam  radius,  and  the 

ratio  of  beam  thermal  energy  to  beam  streaming  energy.  In  this 
case  Rq  » 10  cm  and  the  beam  is  con^sed  of  ions  with  “ 
40  GeV.  16 
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B.  Laser  Absorption  via  Tame  Turbuience 

Recent  experimental  data'  on  absorption  of  laser  light  at  high  intensities  has  confirmed 
the  existence  of  resonant  absorption.  These  experiments  also  indicate  clearly  that  another 
absorption  mechanism  is  present.  For  example,  data  in  ref.  1 shows  that  30%  absorption  occurs 
for  s-polarized  light.  From  this  experiment  and  others,^  we  conclude  that  this  other  absorption 
mechanism  shows  a significant  level  of  absorption,  is  nearly  independent  of  polarization  and 
angle  of  incidence  and  can  co-exist  with  resonant  absorption.  Currently  considered  absorption 
mechanisms  are  listed  in  Table  I along  with  their  characteristics. 

Inverse  bremsstrahlung  falls  off  rapidly  wth  laser  intensity  since  collisions  become  more 
and  more  scarce  as  the  temperature  increases.  Parametric  instabilities  which  occur  near  the 
critical  surface  are  stabilized  by  a short  gradient  scale  length  in  that  region  (due  to  the  deposi- 
tion of  laser  momentum)  and  thus  are  not  expected  to  become  very  important  when  resonant 
absorption  is  significant.  The  most  likely  possibility  then  seems  to  be  ion  acoustic  turbulence 
induced  absorption.  In  addition,  this  mechanism  tends  to  deposit  its  energy  into  electrons  with 
velocities  in  the  range  of  2 to  4 v,„  which  is  also  compatible  with  recent  observations  from 
laser-produced  plasmas. 

Table  I 


Absorption  mechanism 

Absorption 

Angular 

Polarization 

Influence  of  sharp  density 

level 

dependence 

dependence 

gradient  near 

Inverse  bremsstrahlung 

Low 

Weak 

None 

None 

Parametric  instabilities 

? (Probably  low) 

Weak 

Weak 

Unfavorable 

' Ion  acoustic  turbulence 

Moderate  to  High 

Weak 

None 

None 

i 

(Some  in  magnetized  case) 

The  basic  model  used  for  the  study  of  ion  acoustic  turbulence  induced  absorption  is  a set 
of  4 fluid  equations  (density,  velocity,  temperature  and  heat  flux)  coupled  to  laser  light 
(incident  and  reflected)  and  ion  acoustic  modes  equations.  These  equations  have  been  solved 
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in  the  steady  state  starting  from  below  the  critical  density  and  integrating  back  toward  the  laser  * 
Boundary  conditions  have  been  determined  from  the  work  by  Lee  et  al.^  Typical  results 
obtained  from  this  model  are  shown  in  Fig.  Bl.  A high  absorption  level  is  found  which  does 
not  depend  strongly  on  either  polarization  or  target  tilt  angle.  However,  the  main  problem  is 
the  absorption  length  which  reaches  bO/z  in  the  case  of  Fig.  Bl.  It  is  not  at  all  clear  that  the 
plasma  can  expand  that  far  in  case  of  a short  pulse  experiment.  Even  in  a long  pulse  experi- 
ment, a long  underdense  region  could  almost  certainly  give  rise  to  strong  Brillouin  backscatter. 

Although  there  are  many  ways  in  which  this  length  could  be  reduced  (by  assuming  a 
higher  level  of  turbulence  at  x = 0,  by  reducing  the  absorption),  a similar  study  was  under- 
taken in  presence  of  magnetic  held.^  It  is  well  established  that  magnetic  fields  of  order  a mega- 
gauss can  be  generated  in  a laser-produced  plasma. 


In  the  magnetized  case,  the  anomalous  heat  conduction  coefficient  takes  the  form 

\e  <t>  (.k) 


^anom  ^ 


N — r^)  at  X » 0, 


where  k is  the  ion  acoustic  mode  wave  number.  So  for  a given  Q (=  .. 

the  smaller  the  fluctuation  level  is,  the  larger  dTjdx  is  and  therefore  the  larger  the  growth  rate 
is.  Thus,  the  transient  region  of  small  is  effectively  eliminated  in  the  case  of  a magnet- 

ized plasma. 


Also,  the  wave  vector  of  the  ion  acoustic  modes  might  be  aligned  with  the  electric  Held 
vector  and  this  geometry  will  enhance  absorption.  Results  corresponding  to  the  magnetized 
case  are  shown  in  Fig.  B2.  In  that  figure,  20%  absorption  has  been  assumed  at  the  critical  den- 
sity (due  to  resonant  absorption,  for  ex.)  From  these  results,  we  see  that  the  absorption  length 
has  been  reduced  by  a factor  of  more  than  6 and  that  the  absorption  level  is  still  significant. 
The  polarization  and  angular  dependence  are  more  complex  to  determine  than  in  the  unmag- 
netized case  but  are  not  as  strong  as  in  the  case  of  resonant  absorption.  Other  results  show  that 
the  temperaure  ratio  TJT,  can  be  lowered  to  5 in  the  case  of  a magnetized  plasma  for  Z — 2.7 
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(corresponding  to  polyethylene  targets). 


In  summary,  these  results  look  very  promising  and  satisfy  the  characteristics  outlined  in 
Table  I.  One  problem  with  this  approach  however  is  the  boundary  conditions  which  are 
imposed  somewhat  arbitrarily  on  most  fluid  quantities  at  at  • 0 and  on  Q at  jt  - -oo.  The  way 
to  specify  these  quantities  in  a physical  way  is  to  solve  for  the  dynamics  of  the  interaction.  This 
time-dependent  approach  is  currently  being  undertaken  and  should  yield  fully  self-consistent 
boundary  conditions.  The  range  of  integration  needs  to  be  increased  in  order  to  include  both 
the  solid  target  on  one  side  and  the  vacuum  on  the  other  one.  This  model  might  also  shed 
some  light  on  the  time-dependence  of  this  absorption  process. 
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C.  Ablation  Layer  Structure  and  Analysis  of  Rayleigh-Tayior  Instability 

The  NRL  inertial  confinement  program  emphasis  has  moved  toward  pellet  concepts  which 
use  longer  (—  3 to  10  nsec)  lower  intensity  driver  pulses  than  previously  assumed.  For  laser- 
coupled  targets,  this  change  was  motivated  by  our  recent  experiments  with  enhanced  stimulated 
Brillouin  backscatter.  For  electron  and  ion  coupled  targets,  our  motivation  is  the  possibility  that 
10-15  nsec  pulse  lengths  may  soon  be  available.  To  accept  these  10-15  nsec  pulses,  it  is  neces- 
sary to  consider  larger  radius,  thinner  shelled  pellets.  The  computational  studies  of  Rayleigh- 
Tayior  instability  at  NRL  are  being  directed  toward  this  regime. 

There  are  essentially  four  different  circumstances  where  Rayleigh-Tayior  modes  can  grow; 
in  the  ablation  layer,  at  internal  material  interfaces  where  materials  of  different  density  abut,  at 
the  buffer-shell  interfaces  which  occur  in  multiple  shell  pellets,  and  where  the  core  shocks 
interact  with  density  gradients  and  the  pusher  interface.  At  NRL  research  has  proceeded  on 
aspects  of  all  four  of  these  cases.  This  section  of  the  annual  report  is  directed  toward  studies  of 
the  ablation  layer  structure  and  calculations  of  Rayleigh-Tayior  instability  in  the  ablation  layer 
which  have  been  performed  during  the  preceding  year. 

Cl.  Ablation  Layer  Structure 

Previous  numerical  simulations  at  NRL  have  shown  that  density  and  pressure  profiles  for 
an  accelerating  D-T  shell  and  ablation  layer  can  be  temporally  superimposed  with  very  little 
scatter.  Because  of  the  presence  of  these  long  lived  "quasi-static"  profiles,  we  felt  that  it  was 
important  to  model  the  equilibrium  solutions  for  a laser  ablation  layer  both  analytically  and 
numerically  and  to  develop  a detailed  understanding  of  the  nature  of  this  flow.  Facility  with 
such  a model  enables  us  to  calculate  the  quantitative  dependence  of  ihe  s.hell  thickness  and 
acceleration,  the  peak  density,  velocity,  and  temperature  of  the  shell,  the  width  of  the 
Rayleigh-Tayior  unstable  region,  and  the  distance  to  the  critical  surface  on  the  total  plasma 
mass,  the  critical  density,  and  the  absorbed  and  reflected  laser  flux.  This  "quasi-equilibrium",  as 
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determined  by  the  input  parameters  of  the  configuration,  can  be  used  as  an  inexpensive  input 
or  driver  for  other  numerical  studies  and  as  an  equilibrium  solution  for  stability  analyses.  Fol- 
lowing is  the  detailed  description  of  the  laser  ablation  layer  included  in  the  numerical  model. 


For  our  analyses  we  use  the  basic  one-dimensional  steady  state  fluid  equations  for  a 
plasma  slab  transformed  into  an  accelerating  frame  of  reference 

continuity:  (pv)  — 0 (la) 


energy:  (P\  + E\  + q)  - pg\  + la  8(x -x^).  (Ic) 

where  r = r/m„  P — pr.  ^ -I-  pv^.  q — — and  g is  the  accelera- 

tion of  the  slab,  la  is  the  absorbed  laser  flux,  !,  is  the  reflected  laser  flux,  and  x^  is  the  location 
of  the  critical  surface.  As  is  evident,  we  are  treating  the  plasma  as  a fluid  with  a y of  S/3  and  a 
classical  plasma  heat  conductivity  with  a coefficient  variation  of  The  laser  deposition  is 
represented  as  a 6-function  at  the  critical  surface.  Note  that  equation  (la)  tells  us  that 
pv  - constant.  To  facilitate  the  analysis  it  is  convenient  to  transform  the  basic  equations  (1), 
except  at  the  critical  surface,  into  a dimensionless  form  where  the  normalization  variables,  sub- 
script 0,  are  those  defined  at  the  point  of  maximum  density  in  the  plasma  slab. 


The  dimensionless  variables  are 


^({)  - t/t„,  tj({)  - (p„/p)^  £ - gho  X, 
and  the  dimensionless  coefficients  are 


-V..  - 


PoV„ 


km;"^gTl'^ 


and 


The  isothermal  mach  number  is  defined  by 


(3) 


(4) 


(5) 


Only  two  of  the  three  equations  (2)  are  independent.  Depending  on  the  circumstances, 
we  found  that  various  pairs  of  the  three  equations  were  more  convenient  to  analyze  than  oth- 
ers. Linder  this  transformation  the  problem  was  reduced  to  solving  for  two  dependent,  dimen- 
sionless variables  p (the  temperature)  and  r)  (the  density)  based  on  one  independent  spatial 
variable  £ and  two  dimensionless  coefficients  Ng  (primarily  the  thermal  conductivity  usually  a 
large  number),  and  Mg  (the  isothermal  mach  number  at  the  density  peak,  usually  a small 
number).  We  utilized  this  dimensionless  set  of  equations  (2)  to  generate  either  general  numer- 
ical or  analytic  solutions  and  then  matched  these  solutions  to  the  specific  physical  constraints 
and  boundary  conditions  to  determine  the  desired  self-consistent  profiles.  The  numerical  solu- 
tions were  obtained  by  a Runge-Kutta  integration  of  the  non-linear  coupled  set  of  diflerential 
equations  (2)  while  the  analytic  solutions  were  obtained  by  making  appropriate  approximations 
to  the  equations  in  the  various  regions  and  then  solving  the  simplified  set. 

For  the  purpose  of  obtaining  general  analytic  solutions,  the  set  of  equations  (2)  is  con- 
veniently divided  into  three  regions.  The  first  region  is  the  accelerating  shell  (£  < 0).  the 
second  region  is  the  vicinity  of  the  density  peak  (£  = 0),  and  the  third  region  is  the  blowoff 
plasma  (£  >0).  Because  of  certain  basic  physical  trends  for  the  solutions  /3  and  t)  in  each  of 
these  regions,  it  is  possible  to  make  simplifying  assumptions  that  allow  analytic  solutions  to  the 
equations  (2).  As  will  be  discussed  later,  there  are  other  admissible  solutions  to  these 
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equations  with  counterintuitive  physical  trends  but  as  yet  these  have  not  been  sufficiently  inves- 
tigated to  determine  whether  they,  too  are  physically  real. 

a.  The  Region  of  the  Accelerating  Shell  (£  < 0) 

Rewriting  equation  (2a)  into  a more  convenient  form  we  find  that: 


5(^-1)  -I-  - 2£  - - l|- 


(6) 


The  traditional  physical  profiles  for  this  region  require  that  < 0,  0 < 1,  and 

in  words,  the  density  and  temperature  should  decrease  away  from  the  density  maximum  at 
£ — 0 toward  the  inner  regions  of  the  accelerating  shell  (£  < 0),  and  for  this  region  the  flow  is 

isothermally  subsonic.  Therefore  from  (2b)  ^ 1 and  correspondingly 


From  (6)  and  (7)  we  then  obtain  the  following  two  inequalities 

50-1)  -t-M^fr,-!)  -2£  < 0 
and 


(7) 


(8a) 


(8b) 


5(^-1) -t- M.Htj-I) -2£ -1- 

Except  for  a very  narrow  region  where  |£|  ~ is  usually  large),  Eqs.  (8a)  and  (8b) 

imply  that 

5(/3- 1) -1- A/'*(t,-1) -2£  =0  (9) 

is  quite  a good  approximation.  By  differentiating  (2a)  we  find  that 

+ m2  ^ . (10) 

^ rf£  ^ t/£  N„  di  f rf£ 

Again,  except  where  |£|  is  on  the  order  of  -jJ-  and  ^ is  varying  rapidly,  we  have  the 
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approximation  that 


fl{  “C 

Combining  Eqs.  (2b)  and  (11)  gives  us  the  further  approximation  that 

■—  = 0 


1. 


This  latter  result  can  also  be  obtained  from  equation  (2c). 


Collecting  our  results,  we  have  as  a general  approximate  solution  to  the  set  of  equations 


(2)  for  the  region  { < 0, 


5(/3-1)  + A/.2(»,-1)-2{  = 0 


1- 


It  is  interesting  to  note  that  this  approximation  is  independent  of  Ng  and  therefore  the  flow  and 
profiles  in  the  accelerating  shell  are  not  basically  governed  by  thermal  conductivity.  Also  equa- 
tion (13b)  states  that  the  entropy  Pfp^  is  approximately  constant  throughout  the  accelerating 
shell  which  behaves  like  an  adiabatic  gas. 


From  equation  (13b)  we  obtain 


dr,  ^ 3 d0 
di  p*  di' 


Combining  this  result  with  equation  (11)  yields 


M.  ~ 2/3* 

di  5fi*  - 3M} 


It  is  important  to  note  that  equation  (15)  does  not  hold  when  ( - 0 and  - 1 |{  - 0 is  the 

turning  point  for  the  density.  There  “ 0 and  from  (2b)  ■■  l|.  Therefore  it  is 

k-o  ' 


evident  that  varies  from  1 to  approximately — r 

dC  5-3iWj' 


very  quickly  over  a narrow  region. 


point  is  very  helpful  since  it  gives  us  a convenient  position  to  define  as  the  inner  edge  of  the 
accelerating  shell.  This  definition  is  reasonable  since  no  physical  information  about  the  flow 
could  be  transmitted  across  this  point. 

b.  The  Vicinity  of  the  Density  Peak  ({  0) 


which  by  differentiating  equation  (20b)  twice  yields 

N„>  j (22) 

This  is  not  a very  stringent  criteria  for  N^-  Generally  for  physical  problems  of  interest 
>>  5/3.  It  is  also  important  to  note  that  the  region  where  the  Taylor  expansions  (20)  are 
accurate  is  exceedingly  narrow.  Trying  to  utilize  these  results  for  even  a distance  on  the  order 
of  the  thickness  of  the  ablation  layer  does  not  yield  valid  results.  In  fact,  equations  (2)  are 
included  more  for  completeness  than  for  their  practical  importance. 
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Two  very  important  regions  lie  in  the  blowo/f.  The  first  is  the  Rayleigh-Taylor  unstable 
region  where  the  gradients  of  density  and  pressure  have  opposite  sign.  The  second  is  the  vicin- 
ity of  the  critical  surface  where  the  laser  deposition  occurs.  In  this  region  is  increasing 

a{ 

rapidly  from  its  value  of  1 at  { - 0.  Corresponding,  ^ > > 1 except  very  near  the  density 
peak.  It  is  expected  that  this  region  be  isothermally  subsonic  at  least  until  the  critical  surface  is 
reached.  Therefore, 

I - 2{1  <<  5(/3-l)  (23) 

is  a valid  approximation  except  in  a very  narrow  region  about  ( - 0 where  it  is  in  error  by 
about  40%.  Here  retaining  C >s  important  and  th<^  approximation  is  more  properly  written  as 

<<  5(/3-l)-2{  (24) 


The  error  in  (23)  potentially  could  be  as  great  as  20%  in  the  vicinity  of  the  critical  surface  as 
well.  In  practice  (numerical  integrations)  it  never  becomes  that  large.  Under  these  approxima- 
tions equation  (2a)  becomes  either 


Equation  (2Sa)  can  also  be  solved  analytically  to  yield  an  approximation  near  { - 0 that  is  more 
exact  than  (26),  but  little  is  gained  by  including  it  here  due  to  its  complexity.  Though  it  is  not 
crucial  to  use  this  more  exact  approximation  to  define  the  Rayleigh-Taylor  unstable  region,  its 
use  does  improve  the  results  in  a measurable  way.  For  { > 0 it  is  advantageous  to  write  equa- 


tion (2b)  in  terms  of  the  pressure  where 


Po  y/v 


This  gives  two  interesting  forms. 


il'"  f]  ■ I'”  '■''’•I  ■ i i 


2-^  — 2-^  \piP 


Subtracting  (28b)  from  (28a)  yields 


The  second  term  on  the  R.H.S.  of  (29)  is  effectively  a coefficient  times  the  rate  of  change  of 
the  pressure.  Conveniently  it  turns  out  that  the  rate  of  change  of  the  pressure  is  significant 
where  the  coefficient  is  small  and  the  coefficient  is  significant  where  the  rate  of  change  of  the 
pressure  is  small.  Therefore  the  term  is  always  small  and  may  be  neglected. 


We  obtain  as  an  approximation 


tl-jS^exp  [2/3  + A/2(t,-1) -2{]  ^ 

■n- 


The  presence  of  - /3^  in  the  denominator  compensates  for  the  additional  root  at  t)  - j3^  intro- 


31 


‘•'4V 


duced  by  subtracting  (28b)  from  (28a).  Equations  (26)  and  (30)  are  then  an  extremely  good 
approximation  to  the  general  solution  for  { > 0.  It  is  also  relevant  to  note  that  these  analytic 
solutions  with  the  modifications  included  due  to  (2Sa)  will  hold  around  i = 0 and  even  be 

valid  when  { is  slightly  less  than  zero  and  is  varying  rapidly  as  noted  earlier. 

o{ 


From  (28b)  we  obtain  two  interesting  results.  First 


and  second 


Combining  (32)  and  (2b)  yields 


or  from  (5) 


P/P,\  - 0 when 
a{  a{ 


\p/p,\  - 0 When  - 1 

dl  I i dl 


(31) 


(32) 


(33) 


(34) 


Equations  (32)  and  (33)  or  (34)  give  us  respectively  the  rate  of  change  of  the  density  and  the 
temperature  at  the  turning  point  of  the  pressure.  It  is  just  this  offset  of  the  density  and  pres* 
sure  peaks  that  bounds  the  Rayleigh*TayIor  unstable  region,  so-called  because  density  and 
pressure  gradients  are  opposed. 


Combining  (33)  with  equation  (2a)  yields  an  equation  defining  the  relationship  of  the 
variables  at  the  pressure  peak 


_L 

2 ^5  2 


+ M<?(r,-l)-2{ 

^0 


(35) 


The  relevant  question  is  whether  equation  (35)  has  a solution  for  all  values  of  and  N„. 
From  Eq.  (32),  in  the  Rayleigh-Taylor  unstable  region 
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(36) 


0 < i2L  < J_ 
«/{  My 


Therefore 


[M}(r,-\)  - 2{]  < 0 or  - 2{  ^ 0 


there  applying  (37)  to  equation  (35)  we  obtain  the  condition  for  the  existence  of  a pressure 
peak  not  located  at  { - 0, 


Ng  Mgt)  2 


also  in  the  Rayleigh-Taylor  region 


or  from  equation  (27) 


^ [P/Po]  >0  or  > 1 
a{  Po 


Combining  inequalities  (38)  and  (40)  requires  at  the  pressure  peak 

< -y  A/i  5(j3-l)  + . (41) 

It  is  tedious,  but  straightforward  to  show  that  there  can  be  no  regime  where  0 — 1 > 0 and  ine- 


quality (41)  can  hold  when 


N„Mg^  < J >/j. 


This  is  a rather  significant  result  since  it  shows  that  not  all  values  of  Mg  and  Ng  are  admissible 
if  one  requires  a solution  to  the  set  of  equations  (2)  where  there  is  a turning  point  in  the  pres- 
sure profile.  Figure  1 displays  graphically  the  allowable  region  of  Mg  - Ng  space. 

The  only  other  important  modeling  necessary  for  the  region  { > 0 is  near  the  critical  sur- 
face. Here  we  integrate  equations  (1)  across  the  critical  surface  to  obtain  the  jump  conditions 
and  utilize  the  density  shelf  criteria  of  Felber'  and  Lee  at  al}  where 


(43) 


yfT-^  2-v/^ 


(\-My 

- 1 


with  the  subscript  minus  referring  to  the  values  just  interior  to  the  jump  and 

/ t2 


Vc 


Pc 


(pf  is  the  critical  density). 


(44) 


We  force  the  temperature  to  remain  constant  across  the  jump  and  the  heat  flux  to  be  zero  out- 
side the  critical  surface.  We  did  not  believe  that  it  was  important  to  model  the  underdense  gas 
outside  the  critical  surface  since  the  density  is  so  low  there.  Felber'  has  chosen  to  model  this 
region  as  an  adiabatic  gas  and  that  would  probably  be  as  good  a choice  as  any  other.  We  do 
require  that  the  flow  out  of  the  critical  surface  be  supersonic  and  that  gives  us  an  additional 
constraint  on  that 

< Ml  < j (45) 

This  constraint  is  also  included  in  Figure  1.  When  Mg  and  Ng  are  chosen  within  the  bounds 

depicted  in  Figure  1 the  solutions  to  equations  (2)  are  profiles  that  can  be  given  physical  explai- 

nations.  For  values  outside  that  regime  equations  (2)  are  still  solvable,  but  it  is  not  clear 

whether  the  profiles  generated  also  have  a physical  significance  or  are  merely  mathematical 

curiosities.  This  question  is  worthy  of  further  investigation.  There  are  two  possible  approaches 

to  obtain  profile  solutions  to  equations  (2)  for  given  specific  boundary  conditions.  First  the 

coupled  set  of  non-linear  first  order  differential  equations  can  be  integrated  numerically  by  a 

Runge-Kutta  method  and  then  the  boundary  conditions  matched.  Second,  the  general  analytic 

approximations  (13),  (26),  and  (30)  can  be  utilized  instead  of  a numerical  integration  and  then 

matched  to  the  specific  boundary  condition.  Conceptually  there  is  not  much  difference  in  either 

method.  Because  of  the  complexity  of  the  analytic  approximations  both  methods  require  the 

use  of  a computer.  The  rnttior  advantage  of  the  analytic  approximation  approach  is  that,  by 

allowing  an  error  in  the  results  of  a few  percent,  significant  computational  time  can  be  saved. 

Also  the  analytic  approach  often  gives  a better  insight  into  what  is  physically  occurring.  The 

tnree  available  boundary  conditions  not  previously  taken  into  account  are  the  total  plasma  mass, 
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the  absorbed  laser  flux,  and  the  reflected  laser  flux.  Certainly  the  original  slab  mass  must  be 
encompassed  by  the  density  profile.  For  our  model  we  assume  that  all  the  mass  lies  between 
the  interior  isothermal  sonic  point  and  the  critical  surface.  Because  the  densities  outside  this 
region  are  exceedingly  small  compared  to  the  peak  density,  this  approximation  should  have  lit- 
tle effect  on  the  results.  Therefore,  if  m is  the  total  slab  mass  per  unit  area. 


or  from  (4) 


m 


m 


Mo 


(46) 

(47) 


In  either  the  Runge-Kulta  method  or  the  analytic  method  this  integration  (46)  could  be  per- 
formed numerically,  but  for  the  analytic  method,  with  certain  approximations,  the  integration 
can  be  performed  analytically.  From  (46) 


m « 


(48) 


Now  from  (13b)  and  (15)  the  first  integral  in  (48)  can  be  approximated  by 


J»0  1 I 

c,  ^ “ *0^ 


A 03/2  _ 3 Mo 


2 /3 


(49) 


The  dominant  contribution  to  the  second  integral  in  (48)  is  in  the  neighborhood  of  { - 0 where 
T)  ~ /3^  is  a reasonable  approximation  (N.B.  the  approximation  tj  ==  /3^  is  only  of  any  practical 
use  for  this  specific  integration  and  should  not  be  used  elsewhere).  The  second  integration  in 
(48)  can  be  approximated  using  (2Sb)  as 
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Of  course,  in  (SO)  /3,  is,  as  yet,  an  undetermined  parameter. 


At  the  critical  surface  we  handle  the  jump  conditions  and  laser  deposition  much  as  Felber’ 
did  where 


6' 


1 + 2 


Ig-^Q 
P-  V- 


1 + 2 


(51) 


and  we  have  used  the  constancy  of  pv,  and  equation  (5).  It  is  also  true  from  reference  2 that 


CPoT.fi, 


(1-^) 

0 


(52) 


By  use  of  equations  (43),  (48),  (51),  and  (52),  we  are  able  to  determine  the  position  of  the 
critical  surface  for  a given  A/„  and  N„.  This  is  accomplished  by  obtaining  the  value  of  the  first 
integral  in  (48)  either  through  a Runge-Kutta  solution  of  equations  (2)  and  a numerical 
integration  from  { - 0 toward  the  isothermal  sonic  point  or  from  equation  (49).  Then  we 
begin  at  { - 0 and  move  in  the  direction  of  increasing  ( again  using  a Runge-Kutta  solution 
and  numerically  integrating  to  obtain  the  second  integral  in  (48)  or  using  the  analytic  approxi- 
mations and  equation  (50).  We  assume  that  at  each  point  we  have  reached  the  critical  surface. 
Then  from  (48)  we  can  obtain  the  value  t„  such  that  all  the  plasma  mass  would  be  contained 
within  the  density  profile  to  that  point.  With  the  t„  equations  (51)  and  (52)  will  determine  the 
corresponding  p„.  This  value  of  p„  may  then  be  used  to  see  if  equation  (43)  is  satisfied.  When 
(43)  is  satisfied,  we  have  reached  the  critical  surface.  Also  we  have  the  values  of  p,  and 
needed  to  transform  the  general  profiles  into  the  specific  profiles  for  the  given  boundary  condi- 
tions. The  other  values  of  interest  v„  and  g can  be  found  from  (4). 


The  preceeding  discussion  has  masked  a rather  curious  fact.  Even  with  all  the  criteria  we 
have  specified,  we  cannot  determine  a unique  solution  based  on  a single  pair  of  and  N„  but 
rather  have  obtained  a whole  family  of  possible  "quasi-equilibrium”  solutions  based  on  values  of 
M„  and  S„  only  constrained  as  in  Figure  1.  Maybe  we  have  not  specified  all  the  available 
boundary  conditions  but  it  is  difficult  to  think  of  another  two,  only  based  on  the  equilibrium. 


s 


that  we  are  at  liberty  to  specify.  More  likely,  what  the  results  are  telling  us  is  that  there  really 
are  various  possible  ‘quasi-equilibrium”  states  and  just  which  one  would  occur  in  a given  physi- 
cal situation  might  be  determined  by  the  neccessarily  time-dependent  evolution  leading  up  to 
that  equilibrium.  Such  an  effect  might  be  an  initial  shock  which  impacts  an  entropy  to  the 
shell.  Others  may  exist. 

Figure  2 displays  the  profiles  of  density,  temperature,  pressure,  and  isothermal  mach 
number  for  the  ‘quasi-equilibrium”  where  - 4x10"^  /„  - 10.2  TW 

/cm^,  /,  - 10.2^^^'/cm^  m = .773 mg/cm ^ and  = 700.  These  profiles  were 

obtained  by  a Runge-Kutta  integration  of  the  equations  (2),  but  curves  generated  through  the 
analytic  approximations  discussed  earlier  are  virtually  indistinguishable.  These  are  rather  typical 
results  and  while  variation  of  Mo  and  No,  within  the  constraints  of  Figure  1,  does  change  the 
profiles  quantitatively  there  is  not  much  qualitative  variation.  The  density,  temperature,  and 
pressure  are  seen  to  decrease  smoothly  away  from  the  density  peak  toward  the  isothermal  sonic 
point.  For  this  case  the  thickness  of  the  shell  is  about  3S  microns  and  the  peak  density  is  about 
.3g/cm^  The  corresponding  temperatures  and  pressures  at  the  density  peak  are  respectively 
approximately  8 x 10'^  cmVsec^  and  2.5  x 10'^  dyne/cm^  The  acceleration  for  the  shell  is 
about  4 X 10'^  cm/sec^.  It  is  impossible  to  see  in  the  figure,  but  the  offset  of  the  pressure  and 
density  peaks  is  about  .1  micron.  The  peak  pressure  is  only  slightly  greater  than  the  pressure  at 
the  density  peak.  In  the  blowoff  plasma  the  density  drops  sharply  while  due  to  the  thermal  con- 
ductivity the  temperature  increases  rapidly.  It  is  interesting  to  note  that  the  decrease  in  density 


is  nearly  balanced  by  the  increase  in  temperature  such  that  the  pressure  remains  relatively  con- 
stant The  critical  surface  is  approximately  100  microns  from  the  density  peak  and  the  critical 
temperature  is  about  4.5  x lO'^cmVsec^.  The  density  and  isothermal  mach  number  Just  inside 
the  critical  surface  are  respectively  approximately  4.5  x 10^  g/cm^  and  .6.  These  profiles  have 
been  used  as  initial  conditions  for  time  varying  models  and  are  seen  to  remain  relatively 
unchanged  for  long  periods  of  time.  Therefore  they  are  truly  "quasi-equilibrium”  solutions. 


Further,  since  the  density  and  pressure  eradients  in  the  constant  adiabat  shell  are  alianed.  the 


shell  interior  is  Rayleigh-Taylor  stable. 

Development  of  this  "steady-state”  model  is  still  continuing.  The  most  important  question 
yet  to  be  resolved  is  how  to  choose  and  ti„.  What  is  desired  is  to  relate  and  to  some 
physical  properties  that  would  remain  constant  for  the  solutions  while  external  parameters  like 
the  plasma  mass  or  the  laser  intensity  are  slowly  varied.  One  such  possibility  is  the  entropy  of 
the  accelerating  shell.  We  have  utilized  one  of  our  "quasi-equilibrium"  solutions  as  an  initial 
condition  to  a time  dependent  simulation  model  and  followed  the  burning  of  the  shell.  During 
the  evolution,  the  shell  entropy  did  remain  relatively  constant.  We  have  also  examined  the 
constancy  of  the  shell  momentum,  but,  as  yet,  have  not  reached  any  definite  conclusions. 

Our  model  has  been  used  as  an  initializer  to  time  dependent  ablation  models  and  as  a 
pressure  driver  for  a laser-target  compression  simulation.  For  this  latter  case  the  "steady-state" 
model  was  used  to  relate  the  variation  of  the  laser  intensity  to  the  pressure  at  the  surface  of  the 
shell.  Generally  we  used  the  analytic  approximations  for  these  studies,  since  it  entails  a vast 
savings  in  computational  time.  We  also  plan  to  utilize  our  model  as  an  equilibrium  for  a 
piggy-back  analysis  of  Rayleigh-Taylor  stability  in  the  ablation  layer  plasma. 
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NO  PRESSURE  PEAK  BELDV  THE  CURVE 


Figure  Cl. 2 Typical  solution  profiles.  Results  are  for  the  case  where 

= 10“^^,  pg  = 4 X iO“"g/cm  , = 10.2  TO/cm^,  Tj.  - 

10.2  TW/cm^,  m = .773  mg/cm^,  M w .07,  and  r<700. 


Cl.  Analysis  of  Rayleigh-Taylor  Instability 

The  theoretical  treatment  of  Rayleigh-Taylor  modes  in  the  ablation  layer  is  complicated  so 
much  by  the  convection,  conduction  and  acceleration  effects  attendant  in  the  ablation  process 
itself  that  analytic  treatments  are  at  best  approximate.  Numerical  treatments  are  complicated  by 
the  requirement  to  treat  widely  disparate  time  and  space  scales  with  nasty  nonlinearities. 
Therefore  little  quantitative  is  known  about  the  linear  phase  in  laser  pellet  configurations  and 
almost  nothing  qualitative  is  known  about  the  nonlinear  phases. 

Our  picture  of  the  external  (ablation  layer)  Rayleigh-Taylor  instability  was  obtained  as  a 
result  of  a number  of  numerical  simulations  using  the  sliding  rezone  Eulerian  FAST2D''^'^ 
code  in  an  accelerating  frame  of  reference.  Our  current  calculations  have  been  performed  in 
slab  geometry  with  a thermal  conduction  law  and  ideal  hydrodynamics.  These  calcula- 

tions show  that: 

A)  The  initial  linear  growth  for  long  wavelengths  occurs  qualitatively  as  predicted  by  the  sim- 
plest linear  theories. 

B)  There  is  a cutoff  at  short  wavelengths  which  appears  to  force  the  growth  rate  y toward 
zero.  The  cutoff  wavelength  is  roughly  2rrA  where  A is  the  thickness  of  the  ablation 
region. 

C)  The  nonlinear  evolution  observed  is  the  formation  of  a droplet  at  the  location  classically 
occupied  by  a spike.  Bubbles  of  hot  plasma  can  form  in  the  ablator. 

D)  Rayleigh-Taylor  modes  with  k — h (the  shell  thickness)  do  not  immediately  fracture  the 
shell. 

E)  Dynamic  stabilization  by  oscillating  the  ablation  layer  appears  to  stabilize  a wide  band  of 


the  Rayleigh-Taylor  spectrum. 


Figure  C2.1  presents  dispersion  relations  for  the  Rayleigh-Taylor  instability  of  a D-T ict  shell 
determined  from  a series  of  FAST2D  calculations/  ^ Growth  rate  y is  shown  as  a function  of 
wavelength  \ for  two  cases.  On  the  right,  Ko  - S.88  x 10^^  the  usual  plasma  thermal  conduc- 
tion for  D-T,  is  used.  On  the  left  the  thermal  conduction  has  been  increased  artificially  by  a 
factor  of  10.  For  both  cases  the  temperature  at  the  critical  region  was  forced  to  be  280  ev  and 
the  cold  shell  was  initially  20^  thick.  The  size  of  the  plotting  symbols  indicates  an  estimate  of 
the  accuracy  with  which  the  growth  rates  could  be  determined  from  the  numerical  results. 

The  symbols  indicate  which  transverse  mode  in  the  computational  system  is  being  studied. 
Ois  the  first  mode,  the  longest  wavelength  that  the  system  will  support.  The  symbol  x denotes 
the  second  mode,  A the  third  mode,  and  □ the  mode  with  four  wavelengths  across  the  system. 

In  the  standard  case  on  the  right,  three  separate  calculations  were  performed  with  1/it  x 
gridding.  This  was  done  to  get  information  on  both  the  long  and  the  short  wavelength  ends  of 
the  spectrum.  The  cutoff  at  short  wavelength  is  clearly  indicated  in  the  data  but  is  difficult  to 
follow  down  to  zero  growth  because  of  the  timestep  limitations  which  result.  Timesteps,  cou- 
pled with  the  extensive  grid  needed  for  much  finer  resolution  calculations  made  the  cost  of 
exploring  the  convergence  prohibitive.  Nonetheless,  Figure  C2.2  contains  data  from  simulations 
with  varying  spatial  resolution  in  both  directions  showing  apparent  convergence. 

The  system  for  these  calculations  was  initialized  using  the  steady-state,  semi-analytic  solu- 
tions developed  by  Joe  Orens  for  proscribing  the  starting  temperature,  density,  and  velocity 
profiles.*  We  have  verified  that  these  solutions  are  compatible  with  FAST2D  except  for  minor 
discrepancies  when  the  ablation  layer  is  insufficiently  resolved  in  the  simulation  code.  Thus  the 
transient  radial  shocks  driven  by  our  previous  approximate  initialization  are  largely  absent  and 
negligible  growth  is  observed  when  no  perturbation  of  the  equilibrium  is  applied. 

The  initial  perturbation  is  applied  to  the  density  profile  at  its  peak  just  behind  the  ablation 

layer.  The  ten  longest  wavelength  modes  are  excited  with  equal  (small)  amplitudes  and  ran- 
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dotn  phases  corresponding  to  a total  initial  density  perturbation  (RMS)  of  about  .1%  of  the 
shell  thickness.  In  order  to  study  the  mode  amplitude  growths  at  later  times,  a relatively  clean 
diagnostic  was  needed.  Some  characteristic  of  the  slab  had  to  be  Fourier  analyzed  as  a function 
of  y in  order  to  measure  the  different  mode  amplitudes.  We  are  now  using  the  integrated  mass 
through  the  shell  rather  than  the  ablation  layer  location  or  a time  and  space  averaged  vorticity 
as  had  been  done  previously.  This  new  diagnostic  is  insensitive  to  acoustic  waves  in  particular 
which  are  always  observed  rattling  back  and  forth  through  our  accelerating  ablator  shell.  Thus 
the  fluctuations  previously  observed  in  the  Fourier-analyzed  mode  amplitudes  are  greatly 
reduced. 


The  integrated  mass  from  the  inside  of  the  shell  to  the  ablation  layer  is  determined  for 
each  of  the  rows  of  computational  cells,  that  is,  for  each  discrete  value  of  y.  Where  the  laser  is 
burning  through  more  quickly  the  shell  will  be  thinner  and  the  integrated  mass  less.  In  regions 
where  ablator  "droplets"  are  forming  the  shell  is  thicker  and  the  integrated  mass  relatively  large. 
To  determine  the  mode  amplitudes  the  integrated  mass  values  as  a function  of  y were  Fourier 
transformed.  The  growth  rates  were  then  determined  by  plotting  a given  mode  amplitude  as  a 
function  of  time  on  semi-log  paper  and  making  an  estimated  straight  line  fit. 


Figure  C2.3  displays  such  a plot.  After  an  initial  transient  period  an  exponential  growth 
period  develops.  As  can  be  seen,  this  new  diagnostic  allows  for  a clean  measurement  of  the 
growth  rates,  say  for  the  time  interval  from  two  to  four  nsec.  The  growth  rates  determined 
from  this  figure,  along  with  others  based  on  similar  calculations,  is  shown  in  Figure  C2.2.  The 
growth  rates  for  the  first  three  modes  in  each  of  four  different  numerical  simulations  are  plotted 
as  a function  of  the  mode  wavelength.  Error  bars  are  estimated  from  fluctuations  in  the  mode 
amplitude  versus  time  plots.  In  the  first  three  runs  the  x resolution  in  the  vicinity  of  the  abla- 
tion layer  was  1.0^  and  stretched  zones  were  used  in  the  low  density  regions.  The  periodicity 
length  L,  was  taken  to  be  20^  (A  points),  40^  (•points),  and  80/i(Dpoints)  for  three  runs.  In 


each  case  20  equally  spaced  cells  were  used.  The  fourth  run  had  the  number  of  x cells  doubled 


to  120  and  the  x resolution  in  the  shell  correspondingly  halved  to  O.S/i.  The  idea  here  was  to 
see  if  the  0.1^  structure  predicted  for  the  steady  state  must  be  resolved  in  the  computations  to 
achieve  convergence  of  Rayleigh-Taylor  behavior. 

The  shaded  area  in  Figure  C2.2  corresponds  to  the  classical  Rayleigh-Taylor  dispersion 
relation  y - y/icg  where  k - The  width  of  the  bar  is  due  to  the  variation  of  g,  the  system 

acceleration,  as  a function  of  time.  Although  a constant  laser  flux  of  10'^  watts/cm^  was  used, 
the  mass  of  the  shell  decreases  because  of  the  steady  ablation.  Hence  the  acceleration  slowly 
increases. 

At  long  wavelengths  the  runs  become  very  long  because  growth  rates  are  low.  Thus 
appreciable  reduction  of  the  shell  thickness  occurs  and  the  notion  of  a fixed  growth  rate  is 
suspect.  Nevertheless  some  evidence  of  accelerative  stabilization  can  be  argued  from  the  devia- 
tions of  y below  the  'Jkg  limit  shown  clearly  in  the  figures.  This  deviation  is  apparently  quite 
small,  however,  compared  to  what  we  would  expect  based  on  recent  theoretical  predictions 
which  attempt  to  relate  the  Rayleigh-Taylor  growth  rate  to  the  thickness  and  strength  of  the 
unstable  region  where  pressure  and  density  gradients  are  opposed.  An  analysis  of  the  prob- 
lem in  terms  of  vorticity  generation  gives  the  following  equation 

for  the  time  evolution  of  the  vorticity.  A simple  estimate  of  the  growth  rates  is 

= - f\L(^dx  (2) 

0 0 ft* 

Ho 

which  assumes  that  the  m^or  effect  of  unstable  shear  (vorticity  in  the  ablation  layer)  is  to 
change  the  angle  between  Vp  and  VP  rather  than  changing  the  profiles  themselves  when  meas- 
ured locally  normal  to  the  ablation  layer.  This  approximation  is  obviously  better  for  long 
wavelength  variations  as  is  the  neglect  of  convective  terms.  In  Eq.  (2)  d is  the  width  of  the 
unstable  region  from  the  peak  of  p„  to  the  peak  of  Pg.  Z.  is  a length  of  order  A (L  = A/2  for 


the  classical  incompressible  problem)  which  can  absorb  the  geometric  complications.  The  func- 
tion ((x)  is  the  eigenfunction,  unknown  in  general  but  assumed  to  be  constant  when  the 
wavenumber  k is  small  and  the  unstable  region  narrow. 


Growth  rates  calculated  using  this  formula  and  realistic  pjx),  Pjx)  profiles  predict  that 

1 

i low  power-low  acceleration  systems  should  be  more  stable  in  reaching  a given  implosion  velo- 

city than  a high-power  system  because  the  thickness  of  the  Rayleigh-Taylor  unstable  region 
goes  to  zero  as  the  acceleration  decreases.  The  numerical  evidence  for  this  accelerative  stabili- 
zation is  not  all  that  clear  cut  for  long  wavelengths.  There  are  two  possible  explanations.  The 
first  is  that  the  theory  is  wrong  because  the  strong  convection  of  fluid  through  the  ablation  layer 

I 

is  not  properly  accounted  for. 

The  second  explanation,  that  simulation  code  calculations  are  wrong,  can  also  be  argued 
for  the  low  power  ablations  under  discussion.  Ablation  layer  analyses  of  varying  sophistication 
performed  on  homogeneous,  rather  ideal  situations  such  as  treated  in  the  simulations 
above,^  generally  agree  that  0.1/x  resolution  or  better  is  needed  to  resolve  the  details  of  the 
Rayleigh-Taylor  unstable  region  for  the  low  temperature-low  z ablations  of  interest.  The  time 
dependent  2D  and  ID  piggy-back  simulation  models  currently  used  to  calculate  the  Rayleigh- 
Taylor  growth  do  not  resolve  these  scales  effectively  because  the  cost  is  prohibitive.  This  view 
might  also  explain  the  ambiguous  no  growth  results  of  Brueckner,  Jorna,  and  Janda’  as  situa- 
tions where  the  growth  was  much  smaller  than  expected. 

Since  neither  explanation  is  compelling,  the  experimental  guidance  is  lacking,  and  a 

definitive  theory  has  so  far  eluded  us,  the  next  steps  are  to  quantify  the  theory  and  remove  the  j 

I 

numerical  restrictions  which  hamper  time-dependent  computations  with  the  necessary  resolu- 
tion. 

Another  potential  cause  of  instabilities  on  the  surface  of  the  shell  arises  because  of  the 
use  of  multi-shell  targets.  In  order  to  suspend  the  inner  shells  from  the  outer  ones,  support 
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membranes  (webs)  are  being  used.  These  membranes  can  cause  inhomogeneities  in  the  den- 
sity distributions.  The  instabilities  arising  from  these  have  also  been  studied  in  an  idealized 
situation. 

As  a first  step,  we  have  used  FAST2D  in  a slab  geometry  mode.  Our  initial  profiles  are 
the  same  as  was  used  for  our  Rayleigh-Taylor  studies.^  In  order  to  model  the  membrane,  we 
have  introduced  an  osculating,  Gaussian  density  perturbation  consisting  of  the  same  material  as 
the  slab  (here,  D-T). 

All  the  FAST2D  runs  were  performed  on  a 60  cell  by  20  cell  mesh.  The  length  in  the  y- 
direction  was  40/t  corresponding  to  a resolution  of  2/i.  In  the  x-direction  a Im  resolution  was 
used  for  the  slab,  and  the  grid  was  stretched  both  behind  and  in  front  of  the  slab  in  the  x- 
direction.  The  laser  flux  was  10'^  w/cm^  and  the  slab  and  membrane  were  both  composed  of  a 
50-50  Tmixture. 

A uniform  amount  of  mass  for  each  y-cell  was  distributed  between  two  x-cells  along  the 
gaussian.  Because  of  finite  resolution  and  the  grid  stretching  for  large  x,  this  uniform  mem- 
brane appears  as  the  irregular  density  perturbation  in  the  upper  left-hand  corner  of  Figure  C2.4. 
The  density  was  chosen  to  be  0.125g/cc  in  the  unstretched  cells  (Ax  - l/i)  as  compared  to 
0.315g/cc  for  the  maximum  slab  density.  The  five  other  profiles  show  the  density  surface  at 
later  times  as  labelled.  One  sees  that  a large  "T-shaped”  irregularity  develops  in  a relatively 
short  amount  of  time.  A detailed  study  of  the  early  timesteps  shows  that  this  mass  flow  arises 
from  thermal  conduction  effects.  Initially,  the  membrane  is  both  colder  and  denser  than  its 
surroundings.  Hence,  the  membrane  heats  and  also  loses  mass.  The  heating  is  much  more 
rapid,  and  this  leads  to  a pressure  profile  that  forces  mass  in  towards  the  center  in  the  y- 
direction  and  out  away  from  the  slab  in  the  x-direction.  This  forms  the  characteristic  "T-shape”. 

In  Figure  C2.5,  the  time  development  of  two  other  cases  is  displayed.  For  the  top 

sequence,  the  membrane  density  has  been  reduced  by  a factor  of  five.  For  the  lower  one,  the 
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original  density  was  used,  but  the  membrane  has  been  moved  5m  away  from  the  slab.  As  can 
be  seen,  the  prominent  T-shape'  still  develops,  albeit,  more  slowly. 

The  first  three  mode  amplitudes  for  two  of  the  support  web  simulations  are  displayed  in 
the  next  two  figures  (C2.6  and  C2.7).  The  first  is  for  the  standard  case  where  a maximum  den- 
sity of  0. 1 2Sgm/cc  was  employed.  Note  that  the  growth  rate  is  essentially  exponential  and  that 
the  relative  magnitudes  of  the  three  amplitudes  are  very  similar  to  our  Raleigh-Taylor  (R-T> 
simulations  where  an  initial  density  perturbation  is  introduced  instead  of  a membrane.  For 
comparison,  the  approximate  k—\  mode  amplitude  for  our  standard  R-T  case  (see  Figure  C2.3) 
is  represented  by  the  solid  line.  The  second  figure  is  the  same  as  the  first  except  that  it  is  for 
the  reduced  density  (pmem  **  Pmem/S)  simulation.  Here  the  magnitude  of  the  mode  amplitudes 
are  comparable  to  our  R-T  ones. 

We  have  shown  in  our  simulations  that  the  use  of  thin  membranes  (webs)  to  support 

inner  shells  of  multi-shell  pellets  can  lead  to  large  defot.. ns  of  the  ablation  layer  which  are 

certainly  undesirable  for  acheiving  thermonuclear  burn.  Hence,  more  detailed  simulations, 
along  with  experimental  studies,  should  be  undertaken  before  a commitment  to  web-supported 
multi-shell  targets  is  made. 

Dynamic  stabilization  of  the  Rayleigh-Taylor  instability  in  the  ablation  layer  has  been 
demonstrated  for  the  case  where  the  laser  intensity  is  varied  periodically.''^  The  combination 
of  this  phenomenon  with  accelerative  stabilization  (if  it  turns  out  to  be  real)  and  the  natural 
cutoff  at  short  wavelengths  should  permit  successful  implosion  of  shells  with  a dynamic  aspect 
ratio  of  20  or  30  to  1 . We  are  currently  beginning  investigations  of  dynamic  stabilization  of  lay- 
ered ablators  where  thin,  higher-z  layers  temporarily  interrupt  the  heat  flow  to  produce  the  con- 
trolled oscillating  acceleration. 

Figures  C2.8-10  summarize  some  of  our  recent  dynamic  stabilization  calculations  using 
the  FAST2D  computer  code.  The  unstabilized  control  calculation  is  shown  in  the  upper  half  of 


each  figure;  the  dynamically  stabilized  system  is  shown  evolving  in  the  lower  half.  The  shell  is 
composed  of  D-T  treated  as  an  ideal  gas  with  y - S/3.  The  initial  thickness  of  the  shell  was 
20/1  and  the  initial  fluid  (plasma)  profiles  were  chosen  close  to  quasi-static  ID  equilibrium  to 
minimize  initial  transients.  The  outer  surface  of  the  shell  was  rippled  with  a linear  combination 
of  the  10  longest  wavelength  modes.  This  perturbation  amounted  to  a 3.4%  RMS  thickness 
variation  of  the  shell.  The  system  was  driven  by  specifying  the  temperature  of  the  blowoff 
plasma  at  the  critical  density  0.004  gm/cm^  to  be  2S0  ev  in  the  control  case  and  to  oscillate 
about  this  value  at  1.2S  cycles  per  nanosecond  with  an  amplitude  of  12S  ev  in  the  dynamically 
stabilized  case. 

The  small  plot  in  the  lower  left  of  Figure  C2.8  shows  the  spatially  averaged  temperature  at 
the  grid  points  on  the  high  density  side  of  the  average  critical  layer  location  plotted  as  a func- 
tion of  time.  It  also  shows  average  shell  velocities  for  the  two  cases  versus  time.  The  stabiliz- 
ing oscillation  was  started  at  1.37  nsec,  a time  just  before  the  unstabilized  case  began  to  show 
organized  linear  growth.  The  unstabilizated  case  is  shown  as  the  dashed  line.  The  unstabilized 
calculation  was  terminated  at  7.2  nsec,  as  shown  but  the  average  shell  velocity  becomes  a mean- 
ingless quantity  somewhat  earlier.  The  stabilized  calculation  ran  somewhat  longer  and  before 
shell  breakup  achieved  a higher  velocity  by  almost,  a factor  of  two.  The  laser  finally  breaks 
through  in  the  longest  wavelength  mode  supported  by  the  system  (40/t)  in  the  "stabilized”  case 
when  most  of  the  ablator  has  been  burned  away. 

The  upper  part  of  Figure  C2.10  displays  the  Fourier  mode  amplitudes  of  the  mass  integral 
through  the  shell.  The  -I-  data  refer  to  the  mode  of  wavelength  40/t  in  the  stabilized  case. 

is  the  amplitude  of  the  Acth  harmonic  so  KT’  corresponds  to  a sinusoidal  variation  of 
± 0.1/i  when  the  shell  density  is  1 gm/cc.  The  0 and  data  points  refer  to  the  first  and  second 
modes  in  the  unstabilized  case.  Here  the  growth  has  continued  unabated  until  mode  ampli- 
tudes of  AA/  > KT^  are  obtained.  The  corresponding  deformation  of  the  shell  is  shown  at  S.8 
nsec  and  7.2  nsec  in  the  upper  portion  of  Figure  C2.9. 


The  calculation  was  performed  using  l/r  resolution  along  the  shell  and  1/i  resolution 
through  the  shell.  The  runs  required  ~ 2500  timesteps  and  ran  for  about  14  minutes  in  the 
case  of  the  40  x 40  - 1600  grid  point  calculations  shown. 

To  perform  these  dynamic  stabilization  calculations  and  to  resolve  the  ablation  layer  struc- 
ture well  enough  to  study  accelerative  stabilization,  our  ID  plus  piggyback  model  is  being  built 
on  a fuly  implicit  hydrodynamics  module  called  ADINC  which  allows  timesteps  much  longer 
than  the  CFL  condition  would  permit.  Thus  0.1/*  cells  with  sound  speeds  of  order  10^  cm/sec 
can  be  integrated  stably  with  bt  — 10""  sec.  The  model  currently  lacks  an  adequate  treatment 
of  strong  shocks  but  is  otherwise  ideal  for  the  purpose  of  investigating  the  dynamics  and  stabil- 
ity of  quasi-static  ablation  profiles.  A multi-material,  many-layer  shell  is  easily  treated  because 
of  the  Lagrangian  description  embodied  for  the  shell  and  layer  interfaces. 

ADINC  has  been  included  in  a ID  package  that  is  presently  being  used  to  study  our 
quasi-static  ablation  profiles.*  Along  with  ADINC,  this  package  contains  a plasma  thermal  con- 
duction routine,  DIFFUl,  and  an  initializer  that  allows  for  multi-material,  many-layer  shells. 
The  thermal  conduction  for  the  mixture  is  handled  correctly,  and  it’s  charge  dependence  is 
included.'” 

One  of  the  most  important  aspects  of  the  design  of  a successful  laser-powered  thermonu- 
clear energy  source  is  the  ablation  layer.  This  region  must  absorb  and  conduct  the  laser  energy. 
This  must  be  done  efficiently  and  in  a stable  manner  in  order  to  allow  for  maximal  heating  of 
the  fuel.  The  thickness  of  the  ablation  region,  its  shape,  and  its  composition  must  be  deter- 
mined by  simulation  and  experiment. 

As  a first  step  in  simulating  burn  through  of  the  ablation  region,  a one-dimensional, 
multi-material,  plasma  thermal  conduction  code  DIFFUl,  has  been  written.  This  has  been  cou- 
pled with  a one-dimensional  hydrodynamics  code  for  detailed  analyses  of  laser  burn  through  of 


thin  membranes. 


DIFFUl  was  first  tested  by  comparing  our  numerical  results  for  a simple  case  that  has  an 
analytical  solution.  This  was  done  for  the  case  of  constant  density,  initial  temperature  and  ther- 
mal conductivity.  The  right  boundary  temperature  was  specified  according  to  the  analytical 
solution: 


r + Vol 
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(3) 


The  thermal  "shock"  that  is  driven  by  this  temperature  profile  has  a constant  velocity  given  by; 


y 


15  pr„ 


(4) 


A comparison  of  the  numerical  results  generated  by  OIFFUl  and  the  analytical  ones  of  Eqs.  (3) 
and  (4)  are  shown  in  Fig.  C2.ll  along  with  the  parameters  used.  Both  the  "shock'v position  and 
temperature  profiles  are  reproduced  accurately. 


In  order  to  test  the  multi-material  capabilities  of  this  code  we  have  repeated  the  thermal 
"shock"  test,  but  with  the  introduction  of  a 1^  thick  membrane  in  the  plasma.  As  seen  in  Fig- 
ure C2.12  this  has  been  done  for  glass  (SPj),  copper,  and  gold,  all  materials  under  active  con- 
sideration for  use.  We  have  assumed  that  all  materials  are  fully  ionized.  As  expected,  the 
highest  Z material,  gold,  causes  the  greatest  change  from  the  pure  D-T  (no  membrane)  case 
that  is  also  shown. 

Although  no  simple  analytical  solution  now  exists  for  this  problem,  the  numerical  results 
certainly  appear  reasonable.  Since  the  velocity  is  proportional  to  which  goes  like  1/Z^  the 
burn  through  is  slowest  for  gold.  However,  at  large  distances  from  the  slab  all  the  results  seem 
to  display  approximately  the  same  velocity  with  a time  delay  based  on  the  material  and  the  slab 
thickness. 


We  have  shown  in  Figure  C2.12  how  the  burn  through  delay  of  a thin  plasma  slab  is 
dependent  on  its  composition.  A new  series  of  runs  was  performed  on  glass  to  see  how  the 
burn  through  rate  depended  on  slab  thickness,  and  to  study  the  effects  of  the  slab  at  large 


distances. 

In  the  accompanying  Figures  C2.13  both  the  "shock”  position  as  a function  of  time,  and 
temperature  profiles  are  plotted  for  three  cases  as  denoted.  As  expected,  thickening  the  slab 
enhanced  the  retarding  effects  on  the  burn  through.  Also,  for  both  slab  thicknesses,  it  can  be 
seen  that  the  differences  persist  at  large  times  and  distances. 

As  can  be  seen  from  simple  analytical  arguments,  the  velocity  of  the  thermal  "shock" 
should  approach  the  no-slab  value  far  away  from  the  slab.  This  can  be  seen  to  be  verified  by 
comparing  the  slopes  of  the  three  "shock"  position  curves  at  large  times.  However,  neither  the 
"shock"  position  nor  the  temperature  profiles  ever  return  to  their  no-slab  values. 

Preliminary  results  using  the  combined  DIFFUl-ADINC  implicit  package  have  been 
encouraging.  For  example,  one  of  our  quasi-static  profiles^  has  been  shown  to  survive  1,000 
timesteps  (—  lOnsec)  with  very  minor  changes  when  simulated  in  a realistic  time  dependent 
way.  Thus  the  simulation  model  and  the  analysis  are  consistent,  the  model  is  accurate,  and  the 
quasi-static  approximation  seems  very  appropriate. 

A direct  energy  deposition  algorithm  has  been  added  to  the  FAST2D  code.  This  has  been 
done  by  adding  a source  term  into  the  thermal  conductivity  routine  DIFFUT.  Specifically,  a 
time-dependent  two  dimensional  flux  has  been  introduced.  The  x-dependence,  ^-dependence, 
and  time-dependence  are  all  handled  separately. 

Presently,  there  are  several  different  options  available  for  generating  the  time  evolution  of 
the  flux.  These  include  linear,  gaussian,  and  exponential  growths.  Of  course,  any  more  com- 
plicated evolution  that  might  be  needed  to  simulate  realistic  situations  (such  as  a pre-burn 
pulse)  can  also  be  handled.  A sinusoidal  time-dependence  can  be  utilized  to  further  study 
dynamic  stabilization  effects  on  the  Rayleigh-Taylor  instabilities. 

We  have  been  performing  all  of  our  simulations  with  uniform  illumination  in  the  y*- 
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direction.  However,  a more  realistic  Gaussian  profile  is  available.  Also,  the  region  of  illumina- 
tion can  be  limited  to  simulate  the  plug  burn-out  experiments  currently  being  performed  at 
NRL.  Finally,  controlled  perturbations  of  the  laser  intensity  can  simulate  inhomogeneities  that 
exist  in  all  laser  pulses,  and  the  instabilities  arising  from  them  can  be  studied. 

It  is  assumed  that  all  of  the  flux  absorption  occurs  in  a region  around  the  critical  density. 
Peril,  (0.004  gm/cc  for  D-T).  We  distribute  the  flux  linearly  in  the  x-direction  in  the  two  cells 
between  which  p - Peni.  When  he  critical  density  goes  across  a cell  boundary,  a very  steep  rise 
occurs  in  the  flux  of  the  newly  illuminated  cell.  A special  algorithm  has  been  written  to  smooth 
out  this  transition  and  hence  dampen  any  thermal  shocks  that  develop  in  non-physically  due  to 
the  code  resolution. 

The  multi-material  capability  of  FAST2D  has  been  partially  implemented.  Up  to  three 
different  materials  can  be  utilized  simultaneously.  The  equation  of  state  can  handle  both  the 
ion  and  electron  partial  pressures.  The  temperature  has  been  converted  from  the  awkward 
cm^ 

— - unit  to  the  more  standard  e.V.  unit  which  is  appropriate  for  mixtures.  We  have 
sec^ 

developed  a simple  formula  for  calculating  the  thermal  conduction  coefficient  of  a mixture. 
The  derivation  starts  from  the  ion-electron  Coulomb  collisions  and  leads  to  a simple  charge  and 
density  dependence  for  the  thermal  conductivity.  This  has  been  implemented,  along  with  the 
charge  dependence  discussed,  for  example,  in  Spitzer. 

The  mtuor  modification  remaining  involves  allowing  for  the  x-  and  >>-dependence  of  the 
thermal  conductivity  when  solving  the  heat  equation  in  DIFFUT.  Based  on  our  experience 
with  DIFFUl,  a ID  plasma  thermal  conduction  code  where  this  change  has  been  implemented, 
this  modification  will  be  straightforward.  An  example  of  the  multi-material  capability  of 
DIFFUl  can  be  seen  in  Fig.  C2.12. 
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SIMULATION  EXPERIMENT  DISPERSION  RELATIONS 


Figure  C2.1  Dispersion  relations  depicting  growth  rate  y versus  wavelength 
A for  the  external  (ablation  layer)  Rayleigh-Taylor  instability 
as  computed  by  the  FAST2D  con^uter  code.  TVo  different  values 
of  the  thermal  conductivity  are  used. 
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Figure  C2.4  Time  evolution  of  a D-T  slab  irradiated  by  10  w/cm  when 

an  osculating  support  web  is  present.  The  maximum  slab  density 
is  .315  g/cc  and  the  maximum  web  density  is  .125  g/cc. 
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C2.ll  A conparison  of  the  numerical  results  of  DlFFUl  with  the 
analytical  solution  for  the  test  calculation  described  in 
the  text. 
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Figure  C2.12  The  "shock"  position  as  a function  of  time  similar  to 
Fig.  C2.ll,  except  that  a ly  thick  slab  of  different 
materials  has  been  Introduced  (shaded  area) . Hie 
materials  used  and  their  densities  are  as  ledielled. 


D.  Raylelgh'Taylor  Instability  at  Internal  Matciial  Interfaces 


Realistic  numerical  calculations  of  the  Rayleigh-Taylor  instability  at  interior  material  inter- 
faces require  accurate  advancement  of  the  material  interfaces.  This  requirement  favors  the  use 
of  a Lagrangian  formulation  because  Eulerian  techniques  must  either  employ  special  procedures 
to  locate  an  interface  or  use  excessively  fine  gridding  to  resolve  it.  Unfortunately  these  pro- 
cedures are  not  only  cumbersome  and  time  consuming,  but  result  in  a numerical  diffusion 
which  can  grow  large  enough  to  mask  the  physical  motion  of  the  interface. 

The  price  that  is  paid  for  ease  in  locating  an  interface  in  a Lagrangian  treament  is  a dis- 
torted or  skewed  grid  in  all  but  the  simplest  of  flows.  As  shown  in  Figure  Dl,  even  simple 
shear  must  result  in  grid  distortion.  Figure  Dla  illustrates  the  initial  grid  to  either  side  of  a 
shear  layer.  In  the  simplest  case  all  the  zones  on  either  side  of  the  shear  layer  are  moving  with 
a uniform  relative  velocity  with  respect  to  the  zones  on  the  opposite  side.  Figure  Dlb  shows 
the  grid  at  a slightly  later  time  when  the  distortion  has  become  appreciable.  Difference  equa- 
tions for  derivatives  at  interface  vertices  no  longer  involve  the  closest  neighbors  of  a vertex 
since  its  former  neighbors  have  been  shifted  by  the  shear.  An  orthogonal  grid  could  be 
recovered  at  this  time  by  replacing  all  skewed  grid  lines  by  vertical  lines.  However,  reconnect- 
ing any  one  grid  line  results  in  one  triangular  cell  and  one  pentagonal  cell  instead  of  two  quadri- 
lateral ones  as  indicated  by  the  dotted  reconnected  line.  Therefore,  successive  grid  line  recon- 
nections can  result  in  a uniform  quadrilateral  mesh  in  this  case  only  for  a infinitely  long 
(periodic)  mesh.  In  general  this  means  that  grid  reconnections  for  a quadrilateral  mesh  must 
be  global  and  can  result  in  a quadrilateral  mesh  only  for  perodic  grid  anomalies. 

In  contrast,  if  a triangular  grid  is  used,  any  reconnected  grid  line  must  result  in  two  trian- 
gles as  shown  in  Figures  Die  and  Did.  In  other  words  grid  anomalies  can  be  resolved  locally 
without  recourse  to  global  rezoning  and  cannot  result  in  changing  the  character  of  grid  zones. 
Vertices  can  also  be  added  or  deleted  where  required  to  maintain  resolution  with  only  local  grid 
changes.  In  other  words,  it  is  possible  to  dynamically  alter  a triangular  grid  at  a given  location 


with  the  effects  of  the  alteration  being  felt  only  by  the  vertices  near  that  location. 

This  fact  opens  up  the  possibility  of  conserving  physical  variables  locally  even  during  grid 
restructuring.  In  the  hydrodynamics  code  SPLISH  this  goal  is  attained  through  the  use  of  con- 
servation laws  expressed  as  integrals  over  regions  about  each  vertex  and  over  triangles  or  pairs 
of  triangles.  Grid  reconnections,  vertex  additions  or  vertex  deletions  result  in  the  unambigious 
definition  of  new  cell  quantities  through  the  use  of  these  integral  invariants.  Regridding  algo- 
rithms can  therefore  be  formulated  which  are  completely  automated  and  local  and  which 
preserve  conserved  physical  variables. 

During  the  past  year  the  regridding  algorithms  in  SPLISH  have  been  completely  rewritten 
to  accommodate  material  interfaces,  since  the  original  code  was  designed  for  a homogeneous 
fluid.  This  extension  is  not  trivial  since  grid  lines  cannot  be  reconnected  if  they  lie  along  inter- 
faces. To  preserve  a regular  grid,  therefore,  judicious  use  must  be  made  of  vertex  addition  and 
deletion  routines  to  add  or  subtract  vertices  from  the  interface  both  to  preserve  resolution  and 
to  relieve  local  stretched  triangles  bounding  interface  grid  lines.  The  routines  necessary  for  this 
additional  logic  have  been  written  and  used  in  numerical  simulations,  which  will  be  described 
below. 

The  tremendous  flexibility  gained  through  the  use  of  a dynamically  restructuring  triangu- 
lar grid  also  has  its  price.  Global  ordering  was  abandoned  to  allow  reconnections,  but  it  is  pre- 
cisely that  ordering  which  is  necessary  in  the  fast  Poission  solvers  available  for  quadrilateral 
grids.  There  has  been  virtually  no  effort  made  in  the  numerical  community  in  constructing  fast 
solvers  for  an  irregular  triangular  mesh,  and  the  triangle  code  suffers  a corresponding  decrease 
in  computational  speed.  A potentially  greater  problem  results  from  the  loss  of  a one-to-one 
correspondence  between  mesh  cells  and  vertices.  As  has  been  shown  before',  this  results  in  a 
grid  which  can  inhibit  the  physical  flow  if  the  grid  variables  are  not  properly  positioned.  Figure 
D2  illustrates  the  problem  for  the  case  of  a simple  vortex.  A vertex  at  the  center  of  the  vortex 
hu  zero  mean  velocity  as  do  vertices  far  removed  from  the  vortex.  Therefore  there  must  exist 


a ring  of  vertices  which  have  a maximum  velocity  and  are  therefore  moving  in  a shear  flow 
with  respect  to  vertices  radially  displaced  from  it.  In  the  figure  a triangular  cell  is  drawn  for  the 
grid.  In  standard  quadrilateral  difference  schemes  the  divergence  in  the  fluid  is  controlled  by 
adjusting  the  pressure  at  the  center  of  the  quadrilateral.  For  an  incompressible  fluid  the  area  of 
the  quadrilateral  is  therefore  maintained  at  a fixed  value.  In  the  triangular  case  there  are  more 
triangular  cells  then  there  are  vertices.  Using  pressures  centered  in  triangles  results  in  having 

too  many  pressures  which  overdetermines  the  vertex  advancement  required  to  maintain  triangle  j 

1 

area.  In  Figure  D2,  a vertex  having  the  maximum  velocity  is  being  forced  toward  a triangle  | 

I 

side.  Any  attempt  to  maintain  the  triangle  area  must  result  in  lowering  the  velocity  of  that  ver-  ^ 

tex  relative  to  the  other  two  triangle  vertices.  This  clearly  perturbs  the  vortex  motion  which 
these  Lagrangian  vertices  are  supposed  to  model.  Therefore  placement  of  pressures  at  triangles 
centers  in  this  formulation  must  numerically  result  in  altering  the  ph  sical  flow.  This  unfor- 
tunate result  can  easily  be  avoided  by  placing  pressures  at  vertices  nstead  maintaining  a 
’vertex  volume"  as  a measure  of  divergence.  Attempts  to  use  trian^  ul.  gnas  with  improperly 
positioned  physical  variables  in  the  past  have  led  to  the  labelling  of  triangular  grids  as  "stifT  due 
to  exactly  this  problem.  In  Figure  D2,  the  vortex  motion  is  modified  into  a flow  more  closely 
resembling  rigid  body  rotation,  since  in  this  case  only  rigid  body  motion  can  maintain  the  trian- 
gle area  exactly.  The  temptation  to  borrow  from  difference  schemes  formulated  for  quadrila- 
teral grids  remains  very  strong,  however,  to  the  extent  that  triangular  codes  are  presently  being 
used  and  constructed  in  the  DOE  community  which  employ  this  improper  positioning  of  physi- 
cal variables. 

The  fact  that  care  in  defining  the  numerical  algorithms  can  result  in  the  physically  correct 
generation  of  vorticity  is  best  shown  by  a current  simulation  using  the  updated  SPLISH  code. 

As  shown  in  Figure  D3,  the  restructuring  grid  makes  it  possible  to  follow  the  Rayleigh-Taylor 
modes  well  into  the  turbulent  regime. 

The  calculation,  shown  for  two  fluids  of  density  ratio  2:1,  is  patterned  on  the  situation  of 
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an  ablator  (below)  accelerating  a heavier  pusher  (above).  The  free  surface  above  the  pusher 
models  the  gas  region  at  the  center  of  filled  microballoons.  During  implosion  these  interfaces 
are  generally  quite  cold  so  the  thermal  conduction  smoothing  effects  should  be  small.  Convec- 
tive stabilization,  as  it  occurs  in  the  ablation  layer,  is  absent  and  the  density  changes  in  time 
will  be  small  as  well.  Thus  the  classical  Rayleigh-Taylor  problem  is  a good  starting  place  and 
the  turbulent  mixing  theories^  should  be  applicable  here  if  anywhere. 

Figure  D3a  shows  the  initial  grid  of  triangles  with  a small  sinusoidal  perturbation  of  the 
interface  at  r “ 0.  The  initial  velocities  (in  the  accelerating  frame  of  reference)  were  all  taken 
to  be  zero.  Figure  D3b  shows  the  flow  field  as  the  instability  enters  its  nonlinear  phase.  The 
vertex  locations  are  plotted  at  successive  timesteps  to  give  streaks  whose  length  is  proportional 
to  the  local  fluid  velocity.  Shear  flow  is  established  as  the  interface  steepens,  giving  rise  to  vor- 
tices centered  on  the  interface  in  a flow  which  has  been  likened  to  the  onset  of  the  Kelvin- 
Helmholtz  instability.^  In  Figure  D3c  the  two  fluids  are  beginning  to  mix  in  the  Rayleigh-Taylor 
vortices  and  by  Figure  D3d  the  vortices  have  entrained  enough  lighter  fluid  to  form  "bubbles’ 
which  are  about  to  be  completely  enveloped  by  the  heavier  fluid. 

Figures  D3a  and  D3d  illustrate  the  computational  grid  being  used  at  the  timesteps  shown. 
Near  the  end  of  the  calculation  (Fig.  D3d)  most  of  the  grid  carried  by  the  lighter  fluid  has 
flowed  out  from  below  the  downward  jetting  heavier  fluid.  Similarly,  verticies  in  the  heavier 
fluid  have  been  deleted  above  the  upwelling  lighter  fluid.  A reasonable  resolution  has  been 
maintained  by  adding  other  vertices,  particularly  along  the  interface.  Originally  the  interface 
was  resolved  with  10  vertices,  but  46  were  required  by  the  end  of  the  run.  New  regridding  rou- 
tines are  still  needed  to  permit  the  attachment  of  fluid  through  an  intervening  thinned  layer,  as 
shown  in  Fig.  D3d.  The  heavier  fluid  has  not  touched  the  bottom  but  rests  upon  extremely 
thinned  and  skewed  triangle  in  the  lighter  fluid.  The  presence  of  these  near-zero  area  triangles 
has  terminated  the  calculation  due  to  violation  of  the  Courant  condition.  Similar  problems  will 
soon  arise  at  the  free  surface  near  the  upwelling  lighter  fluid  and  at  the  "bubbles'  which  are 
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about  to  be  formed  in  both  vortices. 


Figures  D3c  and  D3d  also  illustrate  a potential  jetting  problem  which  occurs  in  thin  layer 
Rayleigh-Taylor  systems.  In  Figure  D3c  the  motion  of  the  heavy  fluid  into  the  downward  jet 
from  both  sides  creates  a dynamic  pressure  maximum  jetting  fluid  upward  as  well.  This  motion 
is  stable  in  the  sense  that  in  Fig.  D3d  the  jet  is  collapsing,  creating  large  amplitude  waves  at  the 
free  surface.  The  stability  of  this  surface  perturbation  is  sensitive  to  the  thickness  of  the  layer 
of  heavy  fluid  and  may  become  unacceptable  for  thinner  layers.  Alternatively,  the  perturbation 
shown  here  may  become  unstable  at  a later  time  when  the  profiles  or  the  acceleration  have 
changed. 

Larger,  more  finely  resolved  calculations  will  be  performed  to  evaluate  the  rate  of  progres- 
sion of  heavy  fluid  through  the  light  fluid  and  the  quasilinear  effect  of  average  interface 
broadening  on  the  growth  of  longer,  potentially  more  disruptive  modes.  We  will  also  try  to 
identify  flow  regimes  in  which  jetting  may  contaminate  the  fuel. 

Work  was  also  begun  during  the  past  year  in  the  formulation  of  a compressible  Cartesian 
code  using  the  same  integral  invariant  techniques  used  in  SPLISH.  The  most  attractive  method 
involves  the  evolution  of  the  divergences  about  a vertex  through  the  continuity  equation  with 
the  densities  determined  from  the  equation  of  state.  Once  again  vertex  pressures  are  used  to 
advance  vertex  positions  required  to  obtain  those  values  of  vertex  divergence.  The  method  in 
particularly  attractive  because  the  pressure  is  well  behaved  across  material  interfaces,  whereas 
the  densities  are  discontinuous.  Standard  numerical  lore  obtains  densities  from  an  updated  con- 
tinuity equation  and  derives  the  pressures  from  the  equation  of  state.  Our  approach  is  there- 
fore opposed  to  standard  techniques,  but  offers  considerable  incentive  in  terms  of  stability  and 
ease  of  formulation  for  the  triangular  grid.  Because  of  its  radicalness,  a one-dimensional  code, 
the  ADINC  module  discussed  in  Section  C,  has  been  written  to  test  the  concept  before  imple- 
menting it  on  the  triangular  grid. 
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Figure  D2  A single  vortex  test  problem  where  the  use  of  trierngle-centered 
pressures  can  lead  to  nonsense  results.  If  the  triangle  eurea 
is  used  to  determine  compression  locally,  then  the  triemgle 
drawn  will  rapidly  develop  a totally  non-physical  pressure, 
nie  end  result  is  to  lock  the  flow  into  a solid  body  rotation, 
nie  MRL  SPLISH  code  uses  uses  vertex-center  pressures  so  the 
"stiff  grid"  problem  does  not  occur. 
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E.  Instability  and  Symmetry  Issues  in  Colliding  Shells 


Los  Alamos  and  LLL  recently  have  kindled  our  interest  in  velocity  multiplication  by  col- 
liding shells'  ^ ^ The  highest-gain  pellets  seem  to  use  the  technique  to  reduce  the  energy 
requirement  for  ignition  of  large  fuel  reservoirs.  The  impulsive  rebound  of  a low  mass  shell 
from  an  imploding  high  mass  shell  can  at  best  double  the  heavy  shell  velocity  at  the  cost  of  a 
low  energy  efficiency.  To  achieve  this  factor  of  two,  the  onrushing  heavy  shell  can  slow  only 
insignificantly  so  the  amount  of  energy  transferred  to  the  inner  shell  will  be  small.  The  spheri- 
cal case  is  particularly  interesting  because  of  the  geometric  convergence.  The  mass  of  the  outer 
shell  cannot  be  effectively  infinite  even  if  it  is  very  thick.  Thus  the  actual  amount  of  energy 
which  can  be  transferred  is  limited  as  well  as  the  efficiency. 

While  the  efficiency  of  energy  transfer  is  important  to  the  performance  of  high  gain  pel- 
lets, we  are  driven  to  tackle  this  problem  for  other  reasons.  These  solutions  will  help  deter- 
mine the  important  physical  regimes  and  the  general  spherically  symmetric,  one-dimensional 
evolution  of  the  colliding  shell  profiles  for  use  as  dynamic  equilibria  in  Rayleigh-Taylor  subility 
studies.  In  the  nearly  impulsive  collision  of  two  shells,  compression  effects  in  the  shell  material 
is  important  and  will  affect  instability  growth.  Even  in  the  idealized  limit  treated  here,  there  are 
important  trade-offs.  When  buffer  gas  pressure  between  the  shells  is  low  initially,  the  * strike’ 
is  nearly  impulsive.  Classical  arguments  suggest  that  f 'JkgUTdt,  the  number  of  exponentia- 
tions of  the  Rayleigh-Taylor  mode  of  wavenumber  k,  is  minimized  by  such  an  impulsive 
transfer  of  momentum  and  energy.  Symmetrization  by  the  buffer  gas  of  perturbations  on  the 
driver  shell  and/or  isolation  of  the  pusher  shell  from  the  driver  depends  on  a more  gentle 
interaction  however.  The  buffer  gas  should  never  compress  so  much  as  to  let  the  driver 
material  reach  the  pusher  on  a free  fall  line.  Thus  there  certainly  seems  to  be  an  optimum  fill 
pressure.  From  the  dynamics  point  of  view  as  well,  there  seems  to  be  an  important  optimum  to 
be  determined.  As  a function  of  buffer  gas  fill  pressure  and  other  parameters  such  as  pusher 

aspect  ratio  and  material  densities,  what  is  the  optimum  energy  transfer  configuration? 
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To  transfer  the  maximum  energy  the  inner  shell  wants  to  be  thick  enough  and  massive 
enough  to  bring  the  outer  shell  to  rest.  To  transfer  the  maximum  energy  density  the  velocity 
should  double,  hence  the  amount  of  energy  transferred  will  be  f^tr  below  maximum.  At  issue  is 
the  appropriate  quantity  to  optimize  and  the  meaningful  parameter  space  over  which  to  optim- 
ize it.  By  making  the  inner  shell  thin  the  velocity  doubling  can  be  achieved  even  though  the 
shell  is  dense  and  hence  the  energy  density  can  be  maximized.  If  the  velocity  of  a dense 
pusher  is  to  be  doubled,  it  cannot  be  too  thick.  Therefore,  thin  shell  stability  problems  may  be 
expected  to  accompany  configurations  which  simultaneously  try  to  achieve  efficient  energy 
transfer  with  high  energy  density.  With  the  possibility  of  using  the  high  thermal  velocity  buffer 
gas  to  symmetrize  the  implosion,  an  additional  set  of  important  considerations  is  added.  An 
initial  buffer  of  relatively  high  pressure  will  keep  the  two  shells  well  separated  in  transmitting 
the  'strike”  but  will  not  get  as  hot  as  the  initially  lower  pressure  buffer.  Thus  a hotter,  nomi- 
nally better  symmetrizing  buffer  gas  necessarily  allows  the  two  shells  to  get  closer  together. 
Furthermore,  if  the  buffer  gas  gets  hot  enough,  the  possibility  exists  of  more  than  velocity  dou- 
bling by  layering  the  outside  of  the  inner  shell  with  an  ablator.  Thermal  conduction  will  then 
transfer  excess  kinetic  energy  from  the  outer  shell  to  the  inner  shell. 

To  work  on  some  of  these  questions  a very  idealized  model  has  been  developed  and  pro- 
grammed. Figure  El  shows  the  simple,  idealized  model  being  treated  by  the  SHLINT  code. 
The  spherical  volume  of  interest  is  divided  into  five  shells.  The  first  shell  on  the  inside  is  an 
ideal  gas  satisfying  an  adiabatic  compression  law.  This  shell  encompasses  the  origin  and  is  of 
low  enough  mass  that  the  kinetic  energy  is  negligible.  The  third  and  fifth  (outermost)  shell  are 
low  mass  adiabatic  gases  also.  The  density  and  pressure  are  assumed  to  be  constant  throughout 
these  three  gas  regions.  The  values  of  density  and  pressure  are  determined  by  the  volume  of 
the  shell  at  any  instant.  The  fifth  shell  is  bounded  on  the  outside  by  a rigid  impermeable  wall. 
The  second  and  fourth  shells  are  incompressible  fluids  with  separate  but  constant  densities. 
These  shells  are  heavy  compared  to  the  three  gas  shells.  This  system  corresponds  to  the  simple 
mechanics  course  problem  of  a series  of  weights  coupled  by  massless  springs.  Here,  however. 


the  spherical  convergence  adds  an  essential  complication  and  the  linear  problem  of  small  dis- 
placements is  of  no  interest.  Despite  the  complications,  some  analytic  results  are  possible  and  a 
simple  computer  program  was  developed  to  integrate  the  resulting  ordinary  differential  equa- 
tions in  other  situations.  A typical  calculation  involves  specifying  an  inward  velocity  for  the 
outer  shell  and  then  solving  the  coupled  system  of  four  ordinary  differential  equations  for  the 
position  and  velocity  of  the  two  incompressible  shells.  The  system  is  not  dissipative  and  no 
provision  is  made  (currently)  for  internal  energy  in  the  heavy  shells  or  kinetic  energy  in  the 
light  shells. 

This  problem  is  a useful  ,uid  instructive  generalization  of  a problem  originally  treated  by 
Lord  Rayleigh,  the  pressure  developed  in  a liquid  during  the  collapse  of  a spherical  cavity.^ 
Zel'dovich  and  Raizer  discuss  this  problem  from  the  viewpoint  of  similarity  theory  and  sum- 
marize some  of  the  non-ideal  effects  for  the  special  case  of  a bubble  in  water.  ^ Robson  has 
considered  pseudo-spherical  and  cylindrical  implosions  of  liquid  metal  liners  for  the  purpose  of 
adiabatically  compressing  a magnetically  insulated  thermonuclear  plasma.^  Current  interest  in 
the  ICF  program  now  encompasses  multiple  shell  pellets  and  colliding  shell  pellets as  con- 
centration on  exploding  pushers  subsides.  This  added  degree  of  complexity  permits  the  con- 
sideration of  all  sorts  of  new  twists  on  pellet  design. 

There  are  only  four  independent  equations,  two  positions  and  two  velocities,  so  the  sys- 
tem is  fourth  order.  There  are  many  more  parameters,  however.  There  are  five  shell  densities, 
positions  and  energies.  Including  three  gas  constant  parameters  gives  a total  of  eighteen.  This 
is  far  too  big  a parameter  space  to  investigate  — at  least  initially.  To  trim  it  down  some  we  can 
set  the  gas  constants  to  S/3  and  ignore  the  mass  of  the  three  gaseous  shells.  By  assuming  that 
the  fourth  driver  shell  is  thick  when  the  "strike”  on  the  inner  pusher  shell  occurs,  the  two  outer 
radii  and  the  outer  shell  pressure  can  also  be  ignored.  This  leaves  nine  parameters.  To  exam- 
ine the  momentum  and  energy  transfer  from  the  driver  to  the  pusher  through  compression  of 


t the  buffer  gas  we  can  also  ignore  the  initial  energy  of  the  two  inner  shells.  The  remaining 
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seven  parameters  are  the  initial  position,  density,  and  thickness  of  the  second  shell,  the  initial 
position,  density  and  velocity  of  the  fourth  shell,  and  the  initial  pressure  in  the  third  buffer 
shell.  These  seven  parameters  are  reduced  to  four  if  the  driver  shell  inner  edge  position,  velo- 
city, and  density  are  used  to  non-dimensionalize  the  equations.  The  remaining  four  parameters 
I are  the  aspect  ratio,  relative  position  and  density  of  the  pusher  shell  and  the  initial  ratio  of 

energy  densities  in  the  buffer  and  driver  shells  when  the  calculation  is  begun. 

Figure  E2  shows  the  results  of  calculations  initialized  which  the  driver  density,  initial 
driver  inner  radius,  and  initial  driver  inner  surface  velocity  can  be  taken  as  1.0,  1.0,  and  -1.0, 
i respectively,  without  loss  of  generality.  The  driver  is  assumed  to  be  very  thick  during  the  time 

when  the  pusher  is  being  accelerated  appreciably  and  the  final  adiabatic  region  is  given  low 
enough  pressure  that  it  does  not  influence  the  driver  shell  dynamics. 

The  location  rj  of  the  outer  surface  of  the  pusher  was  taken  as  0.1  * r^.  This  choice  is 
also  quite  general  for  small  fill  pressures  since  a different  ratio  of  radii  can  be  made  to 
correspond  to  0.1  by  varying  the  fill  pressure  appropriately  and  scaling  to  the  velocity  the  driver 
would  have  at  r^iO)  - 1.0.  This  approximate  similarity  of  solutions  was  tested  and  the  results 

are  indicated  as  "0"  in  the  figure.  The  — z-  — 0.1  case  was  re-run  for  - 2.0  with  initial  velo- 

A 

city  and  fill  pressure  scaled  to  correspond  at  r4  — 1 .0  with  the  values  used  in  deriving  the  given 
curves. 

The  pusher  density  can  be  similarly  scaled,  at  least  for  relatively  thin  shells.  Velocity 
multiplication  in  planar  geometry  depends  only  in  the  mass  per  unit  area.  To  the  extent  that 
convergence  doesn’t  enter,  the  same  situation  will  hold  for  spherical  coordinates.  The  "xxx" 
points  in  Figure  E2  were  performed  with  " 2,  A//?  — 0.05.  These  values  are  nearly  identi- 
cal to  the  pjj  ” 1.  ^R/R  ” 0.1  case. 

At  low  buffer  pressure  the  transfer  of  energy  is  essentially  instantaneous  on  timescales  of 

I 


wivfvw  for  the  subsequent  implosion.  As  the  initial  buffer  pressure  is  increased,  however,  the 
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interaction  time  increases  because  substantial  momentum  is  transferred  from  the  outer  to  the 
inner  shell  at  greater  separation  distances.  If  the  initial  buffer  pressure  gets  too  large,  the 
pusher  will  implode  prematurely.  For  initial  pressures  in  excess  S x 10~^  of  P4;  y},  the  figure 
shows  poor  energy  transfer  (ijr  < O D even  for  rather  thick  or  dense  pushers  which  do  not 
tend  to  run  away  before  collision. 

The  result  of  these  calculations  is  to  focus  attention  on  the  relatively  sharp  maximum  at 
about  50%  for  modest  aspect  ratio  shells  (10-20:1).  Good  energy  efficiency  is  important, 
clearly,  but  so  also  is  the  original  reason  for  colliding  shell  designs,  velocity  multiplication.  The 
effective  mass  of  the  outer  shell  is  always  reduced  severely  by  spherical  convergence.  There- 
fore, unless  the  inner  shell  is  quite  light,  significant  energy  density  concentration  cannot  be 
achieved.  The  next  step  is  to  set  up  one  of  these  near-optimum  configurations  in  the  FASTID 
code  to  study  the  region  near  the  idealized  optimum  when  more  realistic  equation-of -state  phy- 
sics and  finite  thermal  conduction  is  included. 

Velocity  multiplication  in  spherical  coordinates  has  an  interesting  twist  which  arises  from 
geometrical  convergence.  The  outer  shell  is  freely  falling  in  the  Rayleigh  solution^  before  the 
buffer  gas  pressure  becomes  important.  The  inner  surface  is  accelerating  rapidly  inward  while 
C the  outer  portions  of  the  shell  decelerate  slightly  to  provide  the  momentum  and  energy  flow 

required  to  feed  this  acceleration.  If  the  inner  pusher  shell  is  quite  thin,  the  nearly  impulsive 
transfer  of  energy  will  generally  cause  the  outer  surface  of  the  pusher  shell  to  attain  nearly 
twice  the  velocity  that  the  inner  surface  of  the  driver  shell  had  prior  to  the  "strike’.  Velocity 
multiplication  is  a full  factor  of  two  only  when  the  driver  is  not  appreciably  slowed.  Immedi- 
ately after  the  strike  the  buffer  gas  region  expands  and  cools  rapidly  as  the  pusher  moves  radi- 
i ally  inward  faster  than  the  driven  shell. 


I* 
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At  this  point  spherical  convergence  enters  the  picture  with  a vengeance.  As  soon  as  the 
shells  have  separated  far  enough  for  the  buffer  gas  pressure  to  drop,  the  two  shells  are  falling 
inward  independently.  The  outer  edge  of  the  pusher  starts  decelerating  right  away  to  feed 
energy  and  momentum  to  the  inner  surface.  The  inner  surface  of  the  outer  shell,  on  the  other 
hand,  starts  to  accelerate  right  away  because  of  the  convergence.  Under  a range  of  conditions 
the  driver  catches  up  to  the  pusher  again  before  the  pusher  is  forced  to  turn  around.  This  sub- 
sequent  "restrike”  transfers  energy  from  the  outer  shell  to  the  inner  again,  increasing  the  real 
transfer  efficiency  over  what  is  achieved  on  the  first  strike. 

Figure  E3.  considers  an  idealized  version  of  this  phenomenon  where  the  pusher  remains 
thin  until  the  accelerating  driver  restrikes  it.  The  lower  figure  displays  a series  of  idealized  res- 
trike trajectories  and  the  upper  figure  plots  the  maximum  aspect  ratio  versus  velocity  multiplica- 
tion factor. 

Let  r^it)  be  the  inner  driver  radius  and  r,(t)  the  outer  pusher  radius  as  a function  of 
time.  If  the  shells  collide  at  / - 0 with  the  thin  pusher  initially  at  zero  velocity,  after  the  colli- 
sion the  pusher  is  moving  inward  at  — a and  the  driver  at  — VdiO).  Not  much  energy  is 
transferred  to  a thin  pusher  but  the  driver  velocity  can  be  doubled.  As  long  as  the  pusher  stays 
thin  it  follows  the  trajectory 

r,(r)  - r/0)  - « y/0)  • f. 

If  the  pusher  thickens,  the  outer  surface  actual./  decelerates  but  the  inner  surface  also 
rebounds  sooner  because  of  convergence.  Therefore  recompression  of  the  buffer  and  restrike  is 
more  likely  with  a thin  pusher  than  a thick  one. 

The  initial  shell  collision  is  indicated  by  the  •at  r — 0 in  the  lower  figure.  A number  of 
straight  lines  indicating  thin  pusher  trajectories  have  been  drawn  for  - a VjiO)  with 

a * 1.4,  1.6,  1.8,  and  2.0.  The  smaller  • points  in  the  lower  figure  indicate  restrike  locations 
where  the  driver  catches  up  to  the  pusher  for  each  of  the  pusher  velocities  chosen.  Under  the 
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ideal  circumstances  of  a thick  driver  and  a relatively  light  pusher,  an  analytic  expression  for  the 
driver  triuectory  is  also  possible.  When  the  buffer  gas  is  again  at  sufficiently  low  pressure  after 
the  initial  collision,  the  driver  will  free  fall  with  VjU)  r^U)  - const.  This  driver  equation  of 
motion  can  also  be  integrated  analytically  to  give 


“ r^iO) 


1 - y r7'(0)  r 


2/5 


Since  /•^(O)  - r^fO)  to  define  the  collision,  recompression  of  the  buffer  occurs  when 
r^fr)  - rjir)  where  t = K,,(0)//r^(0)  is  a non-dimensional  measure  of  the  closure  time.  The 
velocity  doubled  pusher  reaches  the  origin  at  r O.S  in  the  absence  of  a restrike.  The  driver 
reaches  the  origin  at  t - 0.4  in  the  absence  of  the  pusher.  Clearly  the  restrike  occurs  before 
T — 0.4  when  the  accelerating  driver  surface  overtakes  the  constant  velocity  or  decelerating 
pusher.  The  upper  figure  is  a plot  of  maximum  aspect  ratio  computed  by  finding  solutions  to 
the  restrike  radius  r^ir)  - r^fr)  numerically.  The  volume  inside  this  restrike  radius,  when  dis- 
tributed uniformly  over  the  surface  of  the  initial  pusher  radius  sphere,  defines  a maximum 
aspect  ratio.  Thicker  pushers  must  rebound  before  restrike  occurs. 


A possible  ramification  of  this  restrike  phenomenon  is  the  use  of  an  ablator  on  the  out- 
side of  the  inner  shell.  The  tendency  of  the  outer  shell  to  "dog  the  tracks’  of  the  inner  shell 
might  possibly  be  used  to  keep  the  buffer  gas  hot  and  the  inner  shell  ablatively  accelerating  due 
to  electron  thermal  conduction.  To  make  this  effective  the  inner  shell  would  have  to  start  out 
relatively  thick  so  convergence  could  be  expected  to  decelerate  the  outer  surface  before  appreci- 
able implosion  has  occurred.  Another  application  of  the  restrike  might  be  to  ensure  that  peak 
compression  of  the  main  fuel  reservoir  in  the  driver  during  restrike  is  timed  to  coincide  with 
ignition  in  the  core  to  optimize  fuel  utilization. 

These  arguments  all  ignore  the  effects  of  compressibility  in  the  heavy  shells  and  the 
effective  loss  of  energy  that  goes  with  it.  Shock  heating  can  be  expected  to  be  a mqjor  contri- 
butor if  the  buffer  gas  is  sufficiently  low  density  that  the  strike  is  impulsive.  If  a softer  strike  is 
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attempted  by  increasing  the  buffer  gas  pressure,  the  inner  shell  tends  to  run  in  an  rebound  off 
the  origin  before  much  energy  has  been  transferred. 

Two  further  applications  of  this  simple  model  are  worthy  of  mention.  To  perform 
dynamic  stabilization  calculations  of  ablators  layered  with  light  and  heavy  material  to  alter  the 
shell  acceleration  periodically,  this  model  will  clearly  prove  useful.  Applications  are  also 
planned  in  the  analysis  of  thin  foils  being  accelerated  by  deposited  laser  and  light  ion  energy. 
In  the  next  section  calculations  on  core  density  and  temperature  profiles  will  be  presented 
which  were  driven  by  the  incompressible  pusher  shell  of  this  idealized  model  as  an  outer  boun- 
dary condition. 
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Figure  E3  Calculations  of  driver-pusher  restrike.  Tlie  lower  figure  shows 
a driver  trajectory  (dashed  line)  intersecting  the  straight  line 
trajectories  of  very  thin  pusher  shells.  The  upper  figure  shows 
the  maximrun  aspect  ratio  (thickness)  the  pusher  could  have  emd 
still  undergo  restrike  as  a function  of  the  velocity  multipli- 
cation ratio. 


F.  RaylcIgh'Taylor  Stability  and  Central  Core  Dynamics 


The  symmetry  of  the  void  closure  during  collapse  of  the  pusher  will  be  critical  in  deter- 
mining the  temperatures  and  compressions  that  can  be  achieved  in  the  central  "ignitor*  region. 
Growth  or  damping  of  Rayleigh-Taylor  instabilities  when  core  shocks  interact  with  the  density 
gradients  formed  in  front  of  the  pusher  or  with  the  pusher  interface  itself  will  play  a large  role 
in  determining  the  degree  of  symmetry  that  can  be  achieved  during  void  closure.  Multiple 
shock  reflection  from  the  central  core  and  the  pusher  interface  potentially  give  several  growth 
opportunities  to  the  Rayleigh-Taylor  instatiblities.  On  the  other  hand,  these  multiple  shocks 
reflections  are  needed  to  achieve  a high  temperature  in  the  central  core  region.  Clearly  we 
must  determine  optimum  ways  of  compressing  the  core  region  to  high  temperatures  while 
avoiding  disasterous  Rayleigh-Taylor  instabilities  which  may  destroy  the  symmetry  of  the  implo- 
sion and  the  resulting  geometric  convergence.  The  timing  of  the  shocks  will  be  determined  by 
a combination  of  All  gas  density,  temperature,  initial  void  radius  and  pusher  velocity  as  a func- 
tion of  time. 

There  are  two  regimes  of  core  density  receiving  about  equal  attention  in  the  search  for  a 
viable  ICF  ignitor.  In  the  first  the  core  density  is  so  low  and  the  temperature  so  cold  that  the 
fuel  mass  and  pusher  are  effectively  imploding  on  a vacuum.  In  the  second  regime  a fluid  of 
sufficient  mass  and  energy  fills  the  core  that  unopposed  implosion  of  the  fuel  is  impossible.  We 
have  diercted  our  attention  to  this  second  regime  in  the  conviction  that  behavior  of  an  ignitor 
with  a fluid  or  gas  fill  will  be  more  reproducible  and  predictable  than  one  based  on  ideal 
vacuum  implosion.  Spall  from  initial  shocks  might  introduce  an  unignorable  amount  of 
material  to  the  void,  for  example.  Another  case  would  be  non-cryogenic,  gas-filled  targets. 

We  hope  to  perform  numerical  simulations  of  the  central  "void"  closure  problem  with  a 
focus  on  stability  and  symmetry  issues.  The  final  calculations  will  be  performed  on  FAST2D 
and  the  new  linear/non-steady  equilibrium  code. 


In  order  to  address  some  o£  these  symmetry  and  stability  problems,  we  have  begun  by 
looking  at  the  dynamic  equilibrium  problem  in  one  dimension.  Preliminary  calculations  are 
being  performed  in  ID  using  the  fully  conservative  FCT  code  FASTI D. 

The  simplest  problem  to  study  first  is  the  spherical  collapse  of  a shock  wave.  Using  a con- 
stant pressure  boundary  condition  on  the  surface  of  a fluid  sphere  which  is  moving  radially 
inward  leads  to  the  formation  of  a spherically  convergent  shock  wave  in  the  gas  filled  core 
region.  Since  the  shock  wave  is  accelerating  due  to  spherical  convergence  as  it  approaches  the 
core,  it  moves  faster  than  a limiting  characteristic  so  that  information  from  the  boundary  condi- 
tion becomes  lost.  Under  such  circumstances  a self-similar  solution  is  established  for  the  case 
of  a constant  core  fill  density.  At  sufficiently  small  radius,  when  the  shock  has  become  strong 
enough  that  the  core  fill  pressure  can  be  neglected  compared  to  the  shock  strength,  a self- 
similar solution  due  to  Guderley  exists.'  This  solution  allows  us  to  check  the  accuracy  of  our 
numerical  solution  in  the  immediate  core  region  at  the  time  of  the  first  shock  collapse. 


The  three  conservation  equations  for  perfect  gas  with  y - 1.4  are  solved  using  the  fully 
conservative  FCT  code  FASTI D.  The  calculations  were  preformed  in  spherical  symmetry  using 
100  zones  on  a variably  spaced  grid.  The  grid  moves  dynamically  to  follow  the  most  important 
region  of  the  flow;  however,  the  calculation  is  still  basically  Eulerian.  The  first  20  zones  con- 
tain a region  in  which  the  grid  spacing  expands  by  a constant  percentage  (~  5%).  The  remain- 

I \p 


ing  80  zones  are  spaced  using  a power  law 


R, 


~ Rt,  + ^R 


i-20 
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with  the  exponent  P chosen 


to  match  the  grid  spacing  at  R^,  the  interface  between  the  two  regions.  The  outer  boundary 
is  allowed  to  move  inward  with  the  boundary  flow  velocity  in  a Lagrangian  fashion  to 
model  the  moving  pusher  interface.  The  boundary  between  the  two  grid  regions  Ri,  is  allowed 
to  move  inward  with  the  flow  as  the  shock  passes  over  it.  This  allows  the  gridding  in  the  cen- 
tral region  to  improve  as  the  shock  moves  inward  without  having  to  pay  the  full  price  of  the 
numerical  CFL  conditions  except  when  the  shock  is  near  the  origin.  The  problem  is  started  at 
an  initial  radius  of  1 cm  with  a zone  size  of  31/1  tbe  origin.  A minimum  resolution  of 
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V is  used  at  the  origin  at  the  time  of  collapse  when  the  outer  boundary  is  at  about  O.S  cm. 


The  outer  boundary  is  held  at  a fixed  pressure  and  alowed  to  move  freely  with  the  fluid  at 
the  boundary.  The  initial  conditions  at  the  boundary  are  set  up  to  satisfy  Rankine-Hugoniot 
shock  with  a pressure  ratio  of  /’i//*!  ~ lOHAf  ~ 30). 


The  self  similar  solution  due  to  Guderley'  predicts  that  the  ratio  of  the  pressure  behind 
the  reflected  shock  to  the  pressure  behind  the  incident  shock  at  a fixed  location  is  a constant 

1/  - /cl^ 


(/’,//’,  •«  26).  In  addition  it  predicts  that  P- 


const  where  P is  the  pressure  behind 


the  shock  at  radius  P and  time  r,  is  the  time  of  collapse  of  the  shock  to  the  origin. 


The  pressure  ratio  is  almost  independent  of  the  final  resolution  at  the  origin.  Varying  this 
resolution  from  4fi  to  I6fi  gives  approximately  the  same  result  (P^/P,  ~ 20)  although  the  peak 
pressure  at  the  origin  varies  by  over  an  order  of  magnitude.  This  is  to  be  expected  since  the 
numerical  problem  is  also  self  similar  and  is  independent  of  grid  size.  The  lower  pressure  ratio 
measured  numei  cally  r<^sults  from  a slight  delay  in  the  reflected  shock  leaving  the  origin  due  to 
the  finite  thickness  of  the  numerical  shock.  This  pressure  ratio  can  never  be  corrected  by  an 
increase  in  resolution  since  the  problem  is  self  similar.  The  delay  can  be  reduced  to  zero  by 
improving  the  resolution  however. 

If  a hard  sphere  with  radius  of  order  the  numerical  shock  thickness  is  placed  at  the  origin, 
it  forces  the  shock  to  be  reflected  sooner.  At  several  radii  the  size  of  this  sphere  is  negligible 
and  the  self-similar  solution  is  again  recovered.  Using  spheres  radius  of  30^  and  50^  with 
minimum  zone  size  of  4fi  both  give  a similar  pressure  ratio  (Pr/P,  ~ 26).  A radius  of  lO^t 
(dudding  out  two  or  three  zones  in  effect)  gives  a lower  value  of  ~22.  There  thus  appears  to 
be  a threshold  radius,  proportional  to  the  grid  spacing  at  the  origin,  above  which  the  analytic 
solution  can  be  recovered.  This  dudded-out  sphere,  for  our  FCT  calculation,  seems  to  play  the 
role  that  adjustable  Von  Neuman  coefficients  do  in  other  calculations. 
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The  size  of  this  sphere  and  the  gridding  needed  will  be  determined  by  the  influence  of 
thermal  conduction  which  destroys  the  seif  similar  solution  at  the  origin  in  any  case.  A sphere 
radius  of  less  than  the  dimensions  of  the  thermal  boundary  layer  will  have  a negligible  influence 
on  the  Anal  solution.  Determining  the  size  of  this  region  is  the  next  step  in  this  investigation. 

I»  - ty 

Shown  in  Figure  FI  is  a plot  of  the  parameter  P 5 — for  the  inward  and  outward 

R *■ 

. U - 

shocks.  The  inward  value  has  a theoretical  value  of  4.3  x 10  ^ given  by  P — --j — - 

i\ 

— - — p ,a^  (the  limiting  value  for  a strong  shock).  Until  the  shock  reaches  the  central  hard 
y-1 

sphere  the  solution  for  the  converging  shock  is  the  same  for  runs  with  and  without  the  hard 
core  and  lies  very  close  the  theoretical  value  of  .43  x 10"^  For  the  diverging  shocks  there  is  a 
considerable  difference  in  the  two  runs  although  both  reach  an  asymtotic  form  of  the  similarity 
solution.  The  case  without  the  hard  sphere  center  always  has  too  low  a pressure.  This  results 
from  the  shock  getting  too  late  a start  reflecting  from  the  origin.  Since  it  starts  late,  it  will 
always  be  moving  into  a region  of  higher  density  and  pressure  than  it  would  have  with  an  ear- 
lier start  and  therefore  never  catches  up  in  the  similarity  sense  even  though  the  delay  is  very 
small. 

The  case  of  an  actual  pusher  collapse  is  much  more  complicated  than  the  simple  self  simi- 
lar Guderley  solution  we  have  been  discussing.  Unfortunately  the  Guderley  solution  is  a limit- 
ing one  valid  in  at  most  a small  region  of  the  real  problem.  The  solution  rapidly  becomes  con- 
taminated by  the  influence  of  the  boundary  which  destroys  the  self-similar  solution  unless  a 
self-similar  boundary  condition  is  used.  Thus  this  soution  will  tell  us  very  little  about  the 
behavior  of  the  core  region  except  near  the  instant  of  the  collapse  of  the  first  shock  wave.  It 
will  give  us  little  guidance  as  to  the  nature  of  the  solution  after  multiple  reflections  of  the  shock 
from  the  core  and  the  driving  pusher  which  in  this  case  is  the  outer  boundary. 

It  is  the  multiply  reflected  shocks  of  the  nonideal  case  that  are  of  primary  interest  in 
determining  the  core  conditions  for  ignitor  design.  Since  the  pusher  will  move  at  a velocity 
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supersonic  with  respect  to  the  All  gas,  shocks  will  be  driven  in  toward  the  origin.  Upon 

reflection  these  outward  moving  shocks  will  rebound  from  the  density  gradients  in  front  of  this 

pusher  and  from  the  pusher  interface  itself.  Each  of  these  reflections,  when  transmitted  to  the 

pusher,  will  slow  it  down,  transferring  part  of  its  energy  to  this  reflected  shock.  The  number  of 

reflections  and  strengths  of  these  shocks  will  depend  upon  the  velocity  of  the  pusher  interface, 

initial  radius,  fill  density  and  temperature  of  the  interior  core  material.  For  a sufficiently  cold  ! 

All,  the  initial  core  temperature  ceases  to  be  a parameter  and  only  the  fill  density  and  mass  will 

be  important. 

If  the  density  of  the  pusher  is  much  larger  than  the  fill  gas  density  and  the  pressure  in  the 
fill  region  sufficiently  lower  than  pC^  in  the  pusher,  the  pusher  may  be  reasonably  modeled  by 
an  incompressible  fluid.  A more  realistic  driving  boundary  condition  than  the  constant  pressure 
used  for  the  Guderley  problem  is  obtained  by  surrounding  the  compressible  core  fluid  with  an 
imcompressible  spherical  shell  which  can  be  used  to  represent  a pusher  for  the  core  fluid.  For  a 
sufficiently  thick  pusher  with  a large  kinetic  energy  compared  to  the  initial  energy  in  the  core 
fluid,  the  concepts  of  a similarity  motion  of  the  pusher  discussed  earlier  in  term  of  the  Rayleigh 
problem^  are  valid.  In  this  case  the  motion  of  the  pusher  is  characterized  by  the  density  of  the 
pusher  and  its  kinetic  energy.  As  the  pusher  moves  in,  it  transfers  part  of  its  energy  to  the 
core  fluid  which  in  turn  slows  the  pusher  down. 

The  question  is  to  find  some  optimum  transfer  of  energy  from  the  pusher  to  the  core 
region  and  optimize  the  temperature  in  the  core.  This  optimum  configuration  is  then  a prime 
candidate  for  Rayleigh-Taylor  analysis.  Clearly  one  can  either  heat  the  core  region  uniformly  to 
a lower  temperature  or  obtain  a very  hot  central  core  with  a surrounding  denser  and  cooler 
region.  Such  considerations  may  lead  to  a variety  of  interesting  ignitor  concepts  but  our  pri- 
mary interest  is  to  obtain  the  dynamic  equilibrium  profiles  for  such  configurations  to  perturb 
about.  Rayleigh-Tarlor  instability  is  likely  in  as  the  reflected  shocks  stow  down  the  pusher  and 
the  density  gradients  generated  in  front  of  the  pusher. 
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Calculations  are  done  using  the  incompressible  model  for  the  pusher  acting  on  a perfect 
(y  » S/3)  fluid  as  the  All.  The  FAST  ID  code  is  used  to  compute  the  motion  in  the  fluid  core. 
The  problem  is  started  with  a shell  at  a radius  of  1000^  surrounding  a fluid  of  density 
PcAS  ” 10”^  gm/cc  at  a temperature  of  10"^  ev.  This  corresponds  to  a flil  fluid  energy  of 
^ 2 X 10~^ir  joules.  The  pusher  is  given  20^  Joules  of  energy  which  gives  an  initial  inner  edge 

velocity  of  the  incompressible  shell  of  - — 710  x 10^  cm/sec.  The  incompressible  shell 
serves  as  the  boundary  condition  for  the  moving  outer  edge  of  the  fluid  region  which  is  com- 

I puted  numerically  using  the  fully  compressible  one-dimensional  code.  The  velocity  of  the  inner 

i edge  of  the  incompressible  shell  is  computed  self-consistently  by  subtracting  from  the  kinetic 

i ' 

j energy  of  the  pusher  the  energy  transferred  to  the  core  fluid  which  can  be  computed  at  each 

1 ' 

I time  step.  Initially  due  to  geometric  convergence,  the  inner  edge  of  the  incompressible  shell 

accelates  until  it  reaches  a peak  velocity  of  — 1.2  x 10^  cm/sec  at  cycle  1150  and  t — 126.9 
ns.  There  after  the  transfer  of  energy  to  the  core  fluid  takes  place  more  rapidly  than  geometric 

f 

convergence  accelerates  it  and  thus  the  pusher  slows  down. 

In  Figure  F2  we  show  the  temperature  as  a function  of  radius  and  the  shell  locations  for 
several  times  in  the  calculation  The  radius  is  in  microns  and  the  temperature  in  units  of 
cmVsec^.  The  first  plot  at  r — 124  ns  corresponds  to  the  situation  just  before  the  first  spheri- 
cally imploding  shock  reaches  the  center.  The  second  plot  at  r — 125  ns  shows  the  reflected 
shock  before  it  has  reached  the  pusher.  The  curve  at  r 126.8  ns  shows  the  shock  after  it  has 
been  reflected  from  the  steep  density  gradient  that  has  formed  just  inside  the  pusher.  This 
curve  is  well  after  the  Guderley  solution  has  broken  down  due  to  the  influence  of  the  boundary 
at  the  pusher.  A second  spherical  shock  is  launched  by  the  pusher  at  this  time  which  overruns 
the  already  heated  core  and  heats  it  to  a peak  temperature  of  1.8  x 10*  ev  shown  by  the  126.9 
ns  curve.  Further  heating  by  yet  more  reflections  is  shown  by  127.16  ns  when  the  pusher  has 
moved  inward  to  its  furthest  point.  As  can  be  seen,  it  is  the  second  shock  that  is  primarily 
responsible  for  heating  the  core  region  to  very  high  tempratures.  In  Figure  F3  we  show  the 
core  density  profile  at  the  time  t — 127.16  sec.  The  profile  shows  a low  density  region  at  the 
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center  where  the  temperature  is  very  hith  sunounded  by  a mass  of  higher  density  inside  the 
pusher.  It  is  this  central  core  that  will  hopefully  serve  as  the  ignitor  for  a pellet  design. 

The  density  at  the  pusher  inner  edge  has  reached  about  9 g/sec  which  is  above  the  density 
of  the  incompressible  pusher.  This  indicates  a breakdown  of  the  assumption  of  an  incompressi- 
ble pusher.  The  implicit  ADINC  model  with  a proper  equation  of  state  for  solid  or  liquid  den- 
sities would  provide  a more  accurate  pusher  for  the  central  region.  Nevertheless,  a great  deal 
can  be  learned  from  these  simple  models  about  the  behavior  of  the  interior  region.  The  next 
step  is  to  include  a realistic  equation  of  state  and  include  the  compressibility  in  the  pusher 
which  will  effectively  absorb  some  of  the  energy  that  the  pusher  delivers  to  the  core  fluid. 

The  steep  density  gradient  formed  in  front  of  the  pusher  will  suffer  from  Rayleigh-Taylor 
instability  when  it  is  decelerated  by  the  reflected  core  shocks.  By  appropriate  timing  of  the  for- 
mation of  this  gradient  and  the  reflected  shocks  the  damage  caused  by  these  instabilities  may  be 
minimized.  The  ultimate  temperature  achievable  and  the  density  obtainable  may  depend  on 
how  well  these  instabilities  can  be  controlled  and  the  degree  of  symmetry  obtainable  in  this  core 
implosion. 

References  and  Footnotes 

I 

1.  G.  Guderley,  Li(fifarh(forschung  19(9),  pp.  302-312,  20  October  1942. 

2.  See  the  discussion  and  references  in  section  E on  the  self  similar  Rayleigh  problem  for 
the  collapse  of  a spherical  bubble. 
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Figure  F2  Temperature  as  a function  of  radius  for  several  tines  for  void 
closure  driven  by  an  incompressible  pusher  of  large  mass.  The 
second  shock  wave  does  most  of  the  heating. 
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CORE  DENSITY  PROFILE 


Figure  F3  Density  profile  at  the  time  of  minimum  radius  for  the  problem 
shown  in  the  previous  figure.  A hot  low  density  core  has 
formed  surrounded  by  a cooler  high  density  region. 
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During  FY78  our  studies  of  heavy  ion  beam  transport  in  quadrupole  focussing  systems 
have  brought  a number  of  useful  and  encouraging  results  to  the  HIB  fusion  program.  The 
most  important  of  these  is  the  demonstration  of  conditions  where  space-charge  limited  beams 
can  be  propagated  long  distances  with  limited  emittance  growth.  Because  of  the  wide  variety  of 
phenomena  which  can  be  investigated  and  which  are  important  for  accelerator  system  design, 
the  choice  of  research  topics  in  the  future  depends  on  the  directions  taken  in  the  heavy  ion  pro- 
gram. Several  of  these  directions  are  apparent.  First,  the  study  of  solenoidal  systems  is  impor- 
tant because  solenoidal,  rather  than  quadrupole,  focussing  will  be  used  in  the  lower  energy  sec- 
tions of  the  accelerator.  It  is  also  important  because  experiments  are  being  undertaken  that  will 
use  solenoidal  geometry  to  study  the  fundamentals  of  focussed  transport.  There  experiments 
will  need  guidance  and  theoretical  interpretation. 

Second  in  importance  is  the  study  of  axial  heavy  ion  beam  transport.  Axial  beam  dynam- 
ics, as  in  bunching,  introduces  another  potential  set  of  instabilities  which  can  degrade  beam 
quality,  transport,  and  focussing.  Therefore  our  simulation  efforts  will  concentrate  on  the  two 
dimensional  r-z  system  where  the  shorter  transverse  time  scales  present  the  m«(jor  numerical 
obstacle. 


We  also  plan  to  continue  our  efforts  on  the  related  light  ion  channeling,  bunching,  and 
focussing.  The  advantages  of  light  ion  beams  for  their  use  as  drivers  in  inertial  confinement 
fusion  are  their  ease  and  efficiency  of  production  and  focussing  and  their  efficient  coupling  to 
matter.  The  m^jor  problems  which  remain  to  be  solved  lie  in  between,  namely  transport  and 
bunching  of  the  beam  from  source  to  target  chamber.  It  is  planned  by  the  NRL  ICF  program 
to  focus  the  ions  first  in  a relatively  short  focussing  chamber  and  then  to  transport  them  in  a 
long  plasma  channel.  We  plan  to  investigate  the  physics  occurring  in  the  plasma  channel  as  a 
result  of  the  beam  passage. 
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The  absorption  of  high  intensity  laser  light  by  matter  remains  a m^or  emphasis  of  the 
DOE  ICF  research  program.  Absorption  by  retum-current-driven  ion  acoustic  fluctuations  has 
been  shown  to  be  an  efficient  absorption  mechanism  which  bridges  the  gap  between  experimen- 
tal absorption  measurements  and  resonant  absorption  theoretical  predictions.  Recent  inclusion 
of  the  magnetic  field  effects  on  the  ion  acoustic  absorption  of  laser  light  have  brought  encourag- 
ing quantitative  success  to  our  absorption  predictions.  Therefore,  we  plan  to  pursue  the  study 
of  laser  absorption  by  this  mechanism  somewhat  further  by  beginning  the  work  on  the  time 
dependent  problem.  We  feel  it  is  important  to  resolve  the  ambiguity  encountered  in  the  steady 
state  solution:  the  fact  that  one  of  the  boundaries  had  to  be  chosen  near  the  critical  density 
where  the  variables  are  not  necessarily  well  known.  Further,  some  of  the  complications  exist- 
ing in  the  steady-state  solution  are  removed  from  the  time-dependent  solution  so  the  work 
should  also  help  to  explain  similar  ambiguities  encountered  in  our  ablation  layer  modelling. 

Our  studies  of  Rayleigh-Taylor  instability  in  laser  ablation  layers  and  associated  physics 
have  constituted  the  single  largest  component  of  our  ICF  theory  effort.  In  the  past  our  theoret- 
ical analyses,  code  development,  and  numerical  simulations  have  been  directed  toward  under- 
standing and  predicting  the  onset,  growth  and  nonlinear  "saturation"  of  these  modes.  This 
emphasis  should  continue  into  the  future  and  will  extend  to  related  ablation  layers  driven  by 
direct  ion  deposition.  Experiments  on  foils  driven  by  lasers,  electron  and  ion  beams 
throughout  the  ICF  community  are  beginning  to  produce  some  experimental  data  on  Rayleigh- 
Taylor  growth  in  real  situations.  We  plan  to  be  in  a position  to  analyze  and  interpret  this  data 
by  a combination  of  one-dimensional  and  two-dimensional  simulations  and  ad  hoc  theoretical 
analyses. 

The  main  theoretical  problem  to  be  resolved  is  the  presence  or  absence  of  accelerative  sta- 
bilization in  the  low  power  regime.  Numerical  resolution  in  ablation-layer  R-T  calculations 
throughout  the  community  has  been  too  coarse  to  date  to  resolve  the  ablation  layer  structure. 
Theory,  on  the  other  hand,  has  proven  up  to  suggesting  a mechanism  for  greatly  reduced  R-T 


growth  but  has  been  helpless  to  prove  the  effect.  We  are  marrying  a piggyback  linearized  stabil- 
ity module  to  our  implicit  ID  Lagrangian  hydrocode.  In  this  way  we  should  be  able  to  provide 
the  SxlO~V  resolution  needed  in  the  ablation  layer  to  resolve  the  separated  peaks  of  density 
and  pressure  profiles  when  the  incident  power  is  ~10'^  watts/cm^  without  working  with  unac- 
ceptably short  timesteps. 

We  also  plan  to  merge  the  quasi-static,  analytic  ablation  layer  model  (developed  during 
FY77-FY78)  with  our  implicit  hydrocode  directly  in  an  attempt  to  demonstrate  the  feasibility  of 
bypassing  the  expensive  ablation  layer  and  blow-off  calculation  altogether.  We  have  found  that 
the  initial  transients  during  prepulse  and  start-up  pervade  even  the  quasi-static  solutions.  They 
specify  the  effective  adiabat  of  the  shell  and  the  shell  momentum  relative  to  the  ablation  layer. 
Spall,  fluff,  preheat,  and  start-up  shocks  all  affect  these  conditions  and  thus  will  be  the  subject 
of  some  work  during  the  coming  year  in  an  attempt  to  determine  better  "initial”  conditions  for 
fluid  dynamic  simulation. 

All  these  studies  and  capabilities  will  be  applied,  as  appropriate,  to  the  possibility  of  using 
dynamic  stabilization  to  limit  the  growth  of  R-T  modes  in  the  ablation  layer.  The  piggyback 
module,  on  top  of  our  implicit  ID  Lagrangian  model,  should  also  pay  important  dividends  in 
analyzing  current  and  planned  foil  experiments  and  possible  future  applications  include  the 
analysis  of  Rayleigh-Taylor  instabilities  in  realistic  colliding  shell  configurations. 

Our  work  on  the  stability  of  interior  material  interfaces  has  a rather  clearly  mapped  out 
series  of  milestones  which  have  to  be  accomplished.  The  mqjor  effort  must  be  in  the  continued 
development  of  the  method  to  include  at  least  the  physics  which  is  currently  in  FAST2D.  This 
will  be  a long  term  effort  with  very  large  potential  payoffs  for  stability  analysis  of  thin  shells  and 
layered  material. 


lem  of  an  ideal  heavy  fluid  over  a lighter  fluid  where  the  short  wavelength  saturation  results  in 
a slow  ’turbulent"  interpenetration  of  the  two  fluids.  Our  object  would  be  to  test  the  Russian 
theory  by  simulation  in  a simple  idealized  case  before  all  the  complications  of  heat  flow, 
compressibility,  and  realistic  equations  of  state  have  to  be  taken  into  account. 

The  second  problem  is  the  simulation  of  a local  conflguration  which  can  be  expected  to 
occur  in  a number  of  ICF  situations.  The  problem  is  the  introduction  of  a thin,  high  density, 
layer  between  two  otherwise  stably  stratified  layers.  Deformation  of  the  interface  as  a whole 
involves  moving  the  thin  layer  as  a surface  and  consequently  long  wavelengths  should  be  stable. 
At  short  enough  wavelength  the  upper  and  lower  surface  of  the  thin  heat  shield  should  be  able 
to  decouple  sufficiently  to  allow  growth  of  Rayleigh-Taylor  modes.  Thus  a cutoff  is  expected,  a 
wavelength  beyond  which  longer  wavelengths  are  stable.  The  mqjor  application  is  to  —1  gm/cc 
ablator  materials  accelerating  ~'0.2  gm/cc  fuel  through  a thin,  20  gm/cc  heat  shield. 

We  also  plan  continued  work  on  symmetry  and  stability  aspects  of  colliding  shells  and  pel- 
let core  dynamics  but  work  in  these  areas  will  continue  to  receive  low  priority  during  the  next 
year  relative  to  our  efforts  on  the  ablation  layer  Rayleigh-Taylor  modes  and  the  material  inter- 
face problems  mentioned  above.  This  is  true  for  two  reasons.  First,  the  conditions  which  exist 
in  the  collision  of  two  shells  or  the  core  collapse  and  void  closure  depend  in  a causal  way  on 
what  goes  on  before.  Not  only  the  level  of  asymmetry  and  R-T  mode  amplitudes,  but  the 
dynamically  evolving  equilbria  themselves  also  depend  on  the  dynamics  and  stability  of  the 
ablator,  the  driver  and/or  the  pusher  shells.  Second,  our  overall  level  of  effort  does  not  war- 
rent  concentrated  effort  on  all  aspects  of  Rayleigh-Taylor  instability  simultaneously.  Neverthe- 
less some  work  in  these  two  areas  is  projected.  The  first  step  is  to  And  meaningful  one- 
dimensional equilibria.  SHLINT  and  our  implicit  hydrodynamics  code  ADINC,  described 
above,  will  be  used  to  zero  in  on  dynamic  equilibria  in  colliding  shells.  FASTI  D,  our  one- 
dimensicnal  FCT  shock  code  will  be  used  to  assess  core  density  and  temperature  profiles.  Sim- 
ple calculations  of  the  general  type  described  earlier  should  be  sufficiently  well  developed  by  the 
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end  of  the  year  to  allow  piggyback  calculations  of  R-T  modes  in  colliding  shells.  Our  modest 
development  program  for  the  2D  Lagrangian  and  Eulerian  simulation  codes  during  this  year 
should  put  us  in  position  to  perform  a few  useful  simulations  of  void  closure  asymmetry 
(FAST2D)  and  growth  of  perturbations  in  shell  collision  (SPLISH)  shortly  thereafter. 
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